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Combustion in Liquid-Fuel Rocket Motors’ 


D. B. SPALDING 


(Imperial College of Science and Technology, London) 


SUMMARY: The paper provides a simplified picture of the processes occurring 
in liquid bi-propellant rocket notors. Droplet vaporisation, chemical reaction, 
and drag between gas and droplet are considered for a one-dimensional model. 
Quantitative theoretical results are presented for the variation of droplet size, 
droplet velocity, gas temperature and gas velocity within the chamber, 
particular attention being paid to the calculation of L*. The theory is applied 
to the German V2 rocket motor. Practical conclusions from, and extensions 
to, the theory are discussed. 


1. Introduction 


1.1. IMPORTANCE OF THE SUBJECT 


Although a rocket motor is mechanically simple it is not easy to develop. The 
reasons are three-fold: — 


(i) The material of the motor must operate under conditions of stress and 
temperature not far from those causing failure; so a poor design means 
not just poor operation but destruction of the unit. 


(ii) The size of the motors used necessitates large testing plant, elaborate 
safety precautions, and expensive services; so testing is slow. 


(iii) The high pressures and temperatures in the combustion chamber, 
combined with the shortness of the practicable time of a test, make 
internal instrumentation difficult. It is consequently hard to determine 
why a change in performance has occurred and to diagnose the cause 
of failures. 


The development engineer therefore has to reach the specified performance and 
reliability by making continual small changes to his design in the light of guess-work. 
The higher the percentage of guesses which are correct, the more quickly and 
cheaply is the final product obtained. It is thus important that he should have in 
his mind an ideal picture of the processes which occur in rocket combustion 
chambers; this picture must be sufficiently simple for design decisions to be made 
in its light; but not so oversimplified. as compared with reality, that the decisions 
turn out to be wrong ones. 
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1.2. PURPOSE OF THE PAPER 

The purpose of the paper is to supply such an ideal picture. This involves 
suppressing: part of the truth and high-lighting those aspects of it that appear 
significant. The result is a crude caricature. but may nevertheless reveal something 
of the character of the subject. 


A main reason for the crudity of the picture to be presented is that it is 
intended to make quantitative deductions, i.e. to employ the techniques of mathe. 
matics. This is necessary because the engineer must have an idea of how big a 
performance change is likely to result from an alteration of design if he is to decide 
whether it is worth while. To introduce into the ideal model all the elements 
which intuition perceives would rapidly overstrain our mathematical capabilities, 

The materials from which the picture will be constructed are facts obtained in 
laboratory investigations of gas flow and droplet burning. It should become clear 
that the results of such investigations, when properly interpreted, can assist the 
development engineer. 


1.3. LIMITATION TO BI-PROPELLANT ROCKETS 


Discussion will be restricted to rockets into which are injected two separate 
streams, one of fuel and the other of oxidant. These are the ones currently used 
for most long-range missiles. Solid propellants, and liquid mono-propellants, will 
therefore be excluded. 


1.4. CONTENT OF THE PAPER 
A brief description of a conventional rocket motor is given first, followed by a 
summary of the main facts and requirements which the designer must consider. 


Then, by distortion and simplification, the idealised rocket motor is introduced, 
its features being chosen so that mathematical deductions from the laws of droplet 
burning and drag lead to predictions of its operating characteristics. 

The relations embodying these predictions are then used to illuminate the 
rocket designer’s task. 


NOTATION 
a_ velocity of sound in the gaseous combustion products 


A.,/Avit the ratio of chamber area to area of injection orifice assumed 
“running full” 


B transfer number, i.e. the mass ratio of droplet material to 
surrounding gas material which, in an adiabatic constant-pressurt 
mixing process, will form a gaseous mixture in equilibrium with 
the liquid surface 


c average gas specific heat 


c(T,—T.) enthalpy difference of combustion products between the exhaust 
state and the boiling temperature of the liquid 
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cu/k Prandtl number of gaseous combustion products 


d_ thickness of propellant sheet flowing away from jet impingement 
zone 


E chemical loading factor (see equation (11)) 


G total mass flux of propellant per unit cross-sectional area of 
combustion chamber 


k average gas conductivity 


L 4 —combustion chamber volume necessary for complete combustion 


throat area 


we 


¢” (rate of fuel consumption by chemical reaction per unit volume and 
time) x (calorific value of fuel) 


q”’,, Value of q”’ at peak of reaction-rate curve (sec Fig. 5) 
Q heat transfer to liquid necessary to vaporise unit mass of it 
r,r, droplet radius, and droplet radius at injection, respectively 
S dimensionless drag parameter (see equation (10)) 
t time 
T absolute gas temperature 
T, boiling temperature of propellant 
T, temperature of gaseous combustion products 
u gas velocity 
v velocity of droplet 
VY, velocity of droplet at injection 
x axial distance from injector 
x* axial distance at which all droplets have disappeared 
% mass proportion of the propellant which vaporises slowly 
y specific heat ratio of the gaseous combustion products 
¢=dimensionless droplet size (see equation (5)) 
f& average gas viscosity 
€ dimensionless axial distance (see equation (12)) 
&* value of € at which droplets have finally disappeared 
p/p, density ratio of combustion products to liquid propellant 
p, droplet density 


t “reactedness” of the gas (see equation (9)) 
x dimensionless droplet velocity (see equation (7)) 
Xo the value of x at injector plane 


» dimensionless gas velocity (see equation (6)). 
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TABLE I 


THE V2 Rocket Mortor®? 


Dimensions 
Throat diameter (internal) 15°75 in. 
Combustion chamber diameter (internal) 36°3 in. 
Nozzle exit diameter (internal) 29°1 in. 
Combustion chamber volume 21,000 in.’ 
L* 113 in. 
Propellants 
Liquid oxygen flow rate 152:5 lb. mass/sec. 
Fuel (75 per cent. C,H,OH. 25 per cent. HO) flow rate 123-5 Ib. mass/sec, 
Film cooling flow rate 13 per cent of fuel 


Operating Conditions 


Combustion chamber pressure 220 Ib./in.? abs. 
Average injector pressure differential 34 Ib./in? 
Cooling jacket pressure differential 63 Ib. /in.? 
Combustion chamber temperature 4760°F abs. 
Thrust at sea level 56,000 Ib. force 


Heat Transfer 


Average heat flux through chamber and nozzle wall 1-3 B.t.u./in.? sec. 
Coolant temperature rise 63°F 


2. Description of a Typical Large Rocket Motor 
2.1. PHYSICAL FEATURES 


To fix ideas, we consider the production motor of the German V2 missile. 
Figure 1 shows a cut-away view, while Table I gives some numerical data. This 
motor is rather “fatter” than more modern ones, which have proportions some- 
where between those of Fig. 1 and those of Fig. 2, which represents a more advanced 
motor for the V2 missile. 


The following features stand out: — 


(i) The rocket motor consists of an axially-symmetrical chamber opening at 
one end to a nozzle, usually convergent-divergent. 


(ii) In the chamber end, opposite the nozzle, are many small orifices through 
which the fuel and oxidant (ethyl alcohol and liquid oxygen in the V2) 
are injected. 


(iii) The nozzle and chamber walls are cooled by fuel which is constrained to 
flow along the outside by means of a cooling-jacket. Additional cooling 
is provided by a film of liquid fuel which is caused to flow along part of 
the inside wall also, immediately in contact with the gas. The latter 
method is not always used, however. 
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FiGuRE 1. Cut-away view of V2 rocket motor. 


1—Fuel intake 5—Oxygen shower head injector 
2—Lower head chamber (Individual head 6—Fuel injector nozzles 

element with oxygen feeder line connector) 7—Combustion chamber 
3—Upper head chamber 8—Expansion beads 


4—Nipple for oxygen feed line 


2.2. OPERATING CHARACTERISTICS AND REQUIREMENTS 


Once ignition has taken place, the chamber is filled with gases at high 
pressure (say 15 atm. abs.) and high temperature (say, 2,650°K)t. The streams of 
liquid fuel and oxidant, on injection into the chamber, break up into droplets; these 
vaporise and burn; the resulting combustion products flow steadily outwards 
through the exhaust nozzle. 


The thrust of the motor is highest if the gases entering the nozzle are uniform 
in state; i.e. they must be uniform in fuel/oxidant ratio and temperature, and 
completely burned. Leaving aside the gases immediately adjacent the wall, this 
means that the injector must be designed so as to produce a fine and even distribu- 
tion of the propellants (for unevenness at the injector face is never wholly obliterated 
by turbulent mixing in the chamber, and will be partially reproduced at the nozzle 
entry); further, it is necessary that the droplets should all have vaporised completely 
before leaving the chamber. The latter process will occupy much of our attention 
in this paper. 


tThese are the V2 figures. Modern pressures and temperatures are higher. 
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FicureE 2. Tubular design of V2 rocket engine. 


The chamber and nozzle walls take up a temperature intermediate that of the 
coolant and that of the gases, the exact value depending on the heat transfer 
coefficients on the two sides. To reduce the gas-side heat transfer rate, and so the 
metal temperature, the fuel/oxidant ratio adjacent the wall is often made higher 
than the average, either by adjustments to the injector or by film-cooling, as in the 
V2; the resultant loss of thrust has to be accepted. However, this aspect of 
combustion chambers is not dealt with here. 


In addition to maintaining itself intact and providing the maximum thrust. it is 
important that the combustion chamber should be as small as possible. Not only 
does size influence weight, but a motor with an excessive surface area is hard to 
keep cool. It therefore appears that, as so often in combustion engineering, two 
incompatible requirements have to be fulfilled; for to make the chamber very small 
would result in unburned droplets passing into the nozzle, and so in a loss of thrust. 


6 The Aeronautical Quarterly 


23. 


( 
been 
inject 
fuel 2 
domi 
atten 


chos 


Febr 


| 
| | 
OH 
| | | 
cf 
— — 
| 
| 
t 
prc 


COMBUSTION IN ROCKET MOTORS 


INJEC TOR 

ov 

= —u 

= 

° 

= ! 

0 = 


FiGureE 3. Model liquid-fuel rocket motor. 


2.3. DESIGN VARIABLES 

Once the thrust, the chamber pressure, and the nature of the propellants have 
been chosen, the major decisions which the designer must make relate to the 
injector: number of orifices, size of hole, injection velocity, impingement pattern of 
fuel and oxidant, and so on. Heat transfer and manufacturing considerations often 
dominate these decisions, but the designer still has sufficient freedom to pay 
attention to the requirements of combustion. 


How should the injector be designed. and the combustion chamber shape 
chosen, if combustion is to be completed in the smallest possible space? 


3. The Ideal Rocket Motor 
3.1. PHYSICAL DESCRIPTION 


In order to be able to answer such questions, we consider the idealised rocket 
motor shown in Fig. 3. This is supposed to have the following properties : — 


(i) Geometry: The injector face is plane and normal to the chamber axis; the 
combustion chamber is cylindrical; the nozzle is convergent-divergent. 


(ii) Injector: This is so designed that all the droplets produced, whether of 
fuel or oxidant, have the same diameter and the same axial velocity. 
There are very many injection orifices, uniformly distributed over the 
injector face. 


(iii) Propellants: The fuel and oxidant have physical and chemical properties 
so related that all propellant droplets vaporise, and are accelerated, at 
the same rate. 


(iv) Conditions in the Chamber: These are one-dimensional, signifying that, 
at any plane normal to the axis, there is only one value of the gas velocity, 
one droplet size, one gas temperature, and so on. This implies that the 
distances between the droplets are not so small that their “boundary 
layers” appreciably interact, and that there is sufficient turbulent mixing 
in the gas phase for the gas condition to be uniform at a cross section 
except in the immediate vicinity of the droplets. 


The ways in which this model departs from reality are obvious : —real injectors 
produce droplets of many sizes and velocities; real fuels and oxidants vaporise at 
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Ficure 4. Variation of droplet and gas conditions along length of model rocket motor. 


different rates; and so on. Models which are nearer to reality can be constructed, 
and their theory worked out, but the present one will serve as a first step. We may 
expect that deductions drawn about it will be approximately valid for real 
rockets also. 


3.2. PERFORMANCE 


Figure 4 sketches some of the information which we would like to have about 
the happenings in a given ideal rocket; of particular interest is x*, the length of 
chamber at which all the droplets have disappeared. What is this length? And 
how can we keep it small? 


The forms of the curves sketched in Fig. 4 may be deduced qualitatively by 
means of quite simple arguments. Thus it is certain that the droplet radius decreases 
continuously during its fiight through the gas; correspondingly the mass flow rate of 
gas, which is zero at the injector, increases continuously with the axial distance +. 
If the droplets are injected at finite velocity «,, they are at first slowed down by 
friction with the gas; farther downstream, however, the gas may be going faster than 
the droplets which therefore start to accelerate. 


However, qualitative information is not enough. It is our task to determine 
the exact shape of the curves in any given circumstance. 
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3.3. THE THEORETICAL PROBLEM 


We suppose that, in a given case, the propellant properties, the initial droplet 
radius, r,, and the initial droplet velocity, v,, are known; we have to predict the gas 
and droplet conditions downstream of the injector, and so to determine x*, the 
minimum length of combustion chamber which will give the full thrust. Then, by 
examination of various starting conditions, we shall be able to tell by how much x* 
can be decreased if the fuel is pre-heated, whether a short fat chamber is better 
than a long thin one, and so on. 


The four physical laws which determine the solution are : — 


(i) The droplet vaporisation law relating rate of radius reduction to the 
local conditions; 


(ii) the droplet drag law, relating rate of droplet acceleration to the local 
conditions; 


(iii) the chemical reaction law, relating rate of change of gas temperature to 
local conditions; 


(iv) the law of mass conservation, relating the increase in gas quantity to the 
decrease in liquid quantity. 


All these laws are fairly well established in quantitative terms. 


The first three of the laws are concerned with rates of change. It is evident 
therefore that the problem, when mathematically formulated, will involve differential 
equations. Further examination reveals that a set of simultaneous first-order 
equations will be obtained. 


Another feature of the mathematical situation is revealed when we note that, 
because both droplets and gas travel in the same direction, the conditions at any 
plane distant x from the injector are only influenced by what happens upstream of 
that plane, but not by what happens downstream. Moreover, the droplet size and 
velocity, the gas velocity. and the prevailing pressure, are all supposed to be known 
at the injector plane, where x=0. It follows that the problem is of the “propagation” 
or “initial-value” type, and that the integration of the equations can proceed 
straightforwardly and without trial-and-error. 


It will appear that in some cases of practical importance, but not in all, 
analytical solution of the equations is possible; that is to say that simple algebraic 
formulae can be found for the required curves and for x*. Even when numerical 
step-by-step integration must be resorted to, however, no serious difficulty in 
calculating the curves presents itself. 


3.4. OUTLINE OF THE PROCEDURE 
The steps in solving the problem are as follows: — 


(a) Collect together the formulae expressing the four physical laws, i.e. set up 
the differential equations. 
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(b) Simplify the equations as much as possible by substituting new dimension- 
less variables for groups of physical quantities. This practice makes the 
calculations easier to handle and see through, and also ensures that a 
single calculation provides the solution to a wide range of actual cases. 


(c) Solve the equations by whatever techniques are appropriate. This 
merely involves “mathematical workshop practice” and need not 
necessarily interest the engineer. 


(d) Interpret the solutions, which are relations between dimensionless 
quantities, in terms of actual rockets. 


Since the details of the theory are available elsewhere“, we here merely discuss 
the physical laws briefly, indicate the meanings of the important dimensionless 
quantities, and then directly introduce graphs which represent the solutions. There- 
after we shall be free to concentrate on their practical importance. 


4. The Theory of an Ideal Rocket Motor 
4.1. THE PHysicaL Laws 
4.1.1. The Droplet Vaporisation Law 


Many studies have been made, both theoretical and experimental, of the 
burning of single droplets in gas streams. As a result, the rate of vaporisation can 
be calculated rather easily and reliably. whether the droplet is at rest or in motion, 
whether it is fuel in an oxidising atmosphere or oxidant in a fuel atmosphere, and 
sO on. 


For the simplest case, the formula for the rate of decrease of radius is® 


r CP 


where the symbols are defined in the Notation, and it is assumed that the Reynolds 
number of the droplet relative motion is small, that liquid enters the combustion 
chamber at a temperature close to its boiling point, and that the average gas 
properties are suitably evaluated. 


The only difficulty which may arise in evaluating equation (1) concerns the 
transfer number B. However, in most liquid-fuel rockets, it is sufficiently accurate 
to use the simple formula 


c (T»- Tx) 


B= (2) 


where c(7T,—T,)=enthalpy difference of combustion products between the 
exhaust state and the boiling temperature of the liquid. 


The important conclusion to be drawn from this section is therefore that the 
rate of radius decrease is easily calculated when certain thermodynamic and 
transport properties are specified. 
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4.1.2. The Droplet Drag Law 


The drag of spheres moving through fluids at low Reynolds numbers is related 
to the droplet size and fluid properties by Stokes’s Law. Combining this with 
Newton’s Second Law of Motion, connecting the drag force with the droplet 
acceleration, we obtain 

dv 9 log. (1 + B) 


where the term (1/B) log. (1 + B) is introduced to account, partially, for the reduction 
of drag caused by the outward mass transfer. 


Inspection of equation (3) shows that the droplet acceleration is also easily 
related to local conditions in the ideal rocket. A correction should be applied to 
the equation when the Reynolds number is not very small. 


4.1.3. The Chemical Reaction Law 


The nature of the ideal rocket is such that the whole gas phase is of uniform 
fuel/oxidant ratio and enthalpy. Consequently, in accordance with well-known 
ideas, the volumetric rate of chemical reactions may be related to the local gas 
temperature by a curve of the form shown in Fig. 5. where ‘ 


q” =(rate of fuel consumption per unit volume and time) x (calorific value 
of fuel) 


T=absolute gas temperature. 


The curve of Fig. 5 represents a property of the propellants and the chamber 
pressure which can, in principle, be determined experimentally. We need not for 
the moment be concerned with the fact that reaction rate data are currently 
insufficient in both quantity and accuracy. 


The relation between g’” and T can be expressed as an equation and combined 
with that representing an application of the First Law of Thermodynamics to a 
small element of the gas phase. It emerges that the gas temperature varies along 
the rocket as a result of the reaction rate (which tends to increase it), and of the 
vaporisation (which tends to decrease it). This is expressed by another first-order 
differential equation. Details may be found in Ref. 1. 


4.1.4. The Law of Mass Conservation 


It is clear that the total mass of propellant flowing across any section of the 
combustion chamber is a constant for a given case; diminution in the liquid flow 
rate, as a result of vaporisation, is compensated exactly by the increase in gas flow 
rate. Consistent with the specification of our ideal rocket, we therefore have 


pu=G(1 = 5) ‘ (4) 
1 
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Hence the gas velocity is connected to the local droplet size by a very simple 
relation. It is interesting to note that the droplet velocity v does not appear in 
this relation. 


4.2. DIMENSIONLESS QUANTITIES 


Manipulation of the equations, and the representation of the results in a general 
manner, is facilitated by the introduction of dimensionless variables and constants, 
defined as follows :— 


¢ will, of course, be unity at the injector (where x=0) and zero at a distance 
x* from the injector. 


(6) 


Gas Velocity: o= 


Since the combustion chamber is supposed sufficiently long for all the liquid 
phase to disappear, (i.e. r becomes zero), comparison of equations (4) and (6) shows 
that © equals unity at x*, but equals zero at x=0. 
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The dropiet velocity is made dimensionless in a similar way to the gas velocity. 


Xo» the value of x at the injector plane, will prove to be an important design 
parameter. x, is easily deducible once the injector design and the liquid and gas 
properties are known, from 


Pi Aorit (8) 
rm 
Gas Temperature : (9) 


r is usually called the “reactedness” of the gas. It is unity for fully burned 
gas and zero for unburned propellant mixture in the liquid state. 


Drag Parameter : (10) 


(cu/k) is the Prandtl number of the exhaust gases and usually has a value 
around 0:7. B is the transfer number discussed earlier; it often has values between 
5 and 10. 


S is therefore a dimensionless property of the propellants; its magnitude will 
ordinarily be of the order of 0°5. 


Chemical Loading: 


t 


E 


k {c(T,—Tx)+Q} log. (1+ B) (11) 
Cp, 


where g’’m is the value of g’” at the peak of the reaction-rate curve of Fig. 5. We 
note that E is large for small propellant droplet sizes and for low values of the 
chemical reaction rate. 


Dist ; 
istance € G 


Ill 


xk 
l 
€ is the dimensionless variable representing distance along the rocket. It may 


be regarded as x divided by twice the distance travelled by a droplet moving at final 
gas velocity before the droplet has completely vaporised. 


The solution of the equations will yield ¢*, which corresponds to x*, and is 
the value of ¢ at which the droplets have finally disappeared. The significance of 


tE is the same as L in Ref. 1. Other notation is unaltered. 
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é* is that it is easily related to the “characteristic length” L* of the rocket. 
Lx, defined by 


x — combustion chamber volume 

throat area 

is widely used as an indication of the compactness of a rocket motor design. It is 

often thought to be a property of the propellants, but this will appear to be 

incorrect. 


The relation between L* and €* may be shown to be 
(k/cp,) log. (1 + B) 


(14) 


Interim Remarks 


(a) It is clear from equation (14) that, if the value of €* were known, insertion 
of quantities relating to propellant properties, injection velocity and droplet size 
would immediately give the value of L*. 


(b) Even before the differential equations have been solved, inspection reveals 
that the desired property must be related to the fuel properties and operating 
conditions by 


The nature of the function f can only be determined by solving the equations. 


(c) \, is unaffected by the droplet radius r, and, although r, enters the chemical 
loading term E, it appears likely (and will be confirmed later) that in any case E 
does not greatly affect €*. We can therefore conclude that r, influences L* only 
in the way apparent from equation (14), namely, 


It is immediately clear that L* is not a property of the propellants alone, but 
also of the droplet size. Fine atomisation permits the combustion to be completed 
in a smaller volume. 


4.3. RESULTS FOR THE CASE OF Low CHEMICAL LOADING (E=0) 
4.3.1. General 


Chemical reaction rate per unit volume is known to be highest at high 
temperatures and pressures. Since these prevail in rocket motors, we can expect 
that the chemical loading parameter E of equation (11) will often turn out to be a 
small number. The physical consequence of high reaction rate is that the gas phase 
is everywhere in thermodynamic equilibrium (fully-burned; reactedness s=1). A 
mathematical consequence is that the differential equation governing 7 may now be 
omitted; the problem is thereby simplified. 
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This case (Model Ia of Ref. 1) is first considered. Discussion of the influence 
of chemical reaction rate is postponed until Section 4.4. 


43.2. Formulae 


This is the case for which the solution takes the form of simple algebraic 
formulae. It is found: — 


(a) that the droplet size ( is related to the distance € by 


(b) that the droplet velocity \ is related to the droplet size ¢ by 
X=Do+3/S-3)] O+1-OS/S-3). (18) 
(c) that the gas velocity » is related to the droplet size ¢ by 


(d) that the minimum chamber length ¢* is related to the drag parameter S 
and the injection velocity x, by 


AGT 3$/ 10 


These are illustrated graphically in the figures that follow. 


4.3.3. Conditions in the Rocket 


Equations (17), (18) and (19) have been evaluated for various values of S and 
\»» to illustrate the conditions in the model rocket. 


Figure 6 shows a set of curves for S=0-5 and ,, =0°5, which may be regarded 
as typical of actual rockets. The following points are noteworthy : — 


(a) The gas velocity increases linearly at first, and then at a reduced rate. 
This characteristic is important; it ensures that the pressure falls steadily 
along the chamber, and renders the flow at the walls similar to that near 
the stagnation point of an aerofoil; it must be taken into account when 
calculating the convective heat transfer to the chamber walls. 


(b) The droplets at first decelerate slightly, and later accelerate very rapidly. 
However, it is only during the last stages of their life that the droplets are 
greatly influenced by the drag forces. 


(c) &€*, where the droplet size € becomes zero, has the value of 0:26. This 
therefore is the sought-after number which must be inserted in equation (14) 
in order to yield the value of L*. With its aid, the size of chamber needed 
for a given propellant system can be calculated. 
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FiGURE 6. S=0°5, x,=0°5. 


Figures 7, 8, 9 and 10 show similar curves for different values of x, and §. The 
trends are illustrated by considering extremes. 


Figures 7 and 8 show, still for S=0-5, the effect of reducing x, to zero and of 
increasing x, to unity. The first implies having a very low injection velocity; the 
second implies one equal to the final gas velocity. In practice y, is probably nearer to 
Othanto 1. It is seen that increasing \, increases €* also. Low x, gives small €*. 


Figures 9 and 10 show the influence of S by representing respectively the 
solutions for S=0 (zero friction) and S=00 (infinite friction). It is interesting to 
note that in the S=0o (unrealistic) case, £*=0-30 regardless of x. The case of 
S=0 is, however, probably much nearer to reality than is S=0o. 


4.3.4. The €* Relation 


Figure 11 clarifies the way in which the “minimum length” parameter {* 
depends on the drag parameter S and the injection parameter y,. It is seen that 
increasing the injection velocity always increases the necessary length of chamber. 
but that the effect is less at the larger § values. It is necessary to remember that 
S is a propellant property which will always be around 0-5 and is not at the disposal 
of the designer; also x,=1 means that the injection velocity is equal to the final gas 
velocity in the chamber. 
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FIGURE 7. S=0°5, X,=0. 
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1-0 

W (GAS VELOCITY) 
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Ficure 9. S=0, x,=0°5. 
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FiGurE 10. S=%, all x,. 
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FIGURE 11. 


44. RESULTS WHEN THE CHEMICAL LOADING IS NoT SMALL (E + 0) 


The discussion and results of Section 4.3 hold when the chemical reactivity of 
the propellants is so high that their reaction rate is controlled solely by the physical 
process of vaporisation. We now consider what happens when this is not the case; 
in so doing we shall find out quantitatively when chemical reaction is important and 
when it is not. 


The corresponding difference in the mathematical situation is that a non-linear 
differential equation for the temperature distribution has to be solved in addition to 
the previous ones. The solution procedure is somewhat more laborious, being 
numerical rather than analytical. Here, however, we merely consider sample 
results; full details can be found in Refs. 1 and 4 (Rocket Model II). 


It is convenient to fit an algebraic formula to the reaction-rate expression of 
Fig.5. The formula used here is 


Which has the right sort of shape. @”’:, is the maximum value of the reaction rate 
of the gases, i.e. it is the ordinate of the peak of the curve of Fig. 5. Other 
expressions can be used without greatly altering the results. 
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Ficure 12. 


4.5. CONDITIONS IN THE ROCKET 

Figure 12 shows the results of calculations for a Model II rocket with Xo=05 
and S=1. The method of representation is slightly different from before, the 
abscissa being 1—¢. It will be recalled that ¢ is the dimensionless droplet radius, 
so 1—{ is zero at the injector and unity at the completion of combustion. 


Curves for x, the droplet velocity, and 7, the gas temperature, have been plotted 
for values of E between zero and 0:3. The gas velocity » has been omitted for 
clarity; it would, of course, traverse the rectangle diagonally upwards from right 
to left. 
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FicureE 13. The influence of chemical reaction on €*. 


The most notable feature of Fig. 12 is that the relation between x and 1-—¢ 
depends very little on the value of the chemical loading E. The 7r-curves are greatly 
influenced, however, in the sense that 7 falls at any point as E increases. This is 
to be expected: the chemical reaction rate is now not fast enough to keep the 
temperature at its maximum value in the face of the constant addition of cold 
vapour to the mixture. 


It will further be seen that 7 is not unity at 1—C=1, where the droplets have 
disappeared. Consequently, the chamber will have to have a greater length than 
that given by €* in order that the gases should be fully burned before entering the 
nozzle. Strictly, the gases are only fully burned after an infinite length; however 
the inefficiency, 1—7, will have been reduced to 0-368 of its value at €* at an axial 
distance equal to £* + E. 


46. LENGTH OF CHAMBER 
Since € did not appear explicitly in Fig. 12, the effect of E on the length of 
chamber could not be discussed. This is shown in Fig. 13. 


It appears that £* remains constant at first, as E increases from zero; eventually 
itrises somewhat. 


47. FLAME EXTINCTION 


The most striking consequence of the calculations only emerges from Figs. 12 
and 13 by default; namely, combustion is not possible at chemical loadings greater 
than 0-3 (for x,=0°5). 
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More generally, it is found that steady combustion is only possible if th 
quantity E/ x, is less than a certain number, which depends on the shape of the 
reaction-rate expression but not on S, which is 0-6 for the present example, an{ 
which is never less than 0-33. 


This conclusion brings rocket combustion into line with combustion in the gas 
turbine and ram-jet, where an excessive rate of fuel supply can cause extinction of 
the flame. Inspection of the structure of E shows that excessively fine atomisa. 
tion, or too low an injection velocity, might cause E/ x, to rise above the critical 
value: extinction would ensue, the propellant spray itself having acted as the 
quenching agent. 


It is unlikely that any practical rockets operate with values of E/x, anywhere 
near the critical under normal conditions, and this point is examined later. 


5. Discussion and Extension 


The quantitative implications of the foregoing results for the rocket designer 
are now considered. This can be done by considering data typical of actual rocket 
motors. In some respects, actual motors differ considerably from the ideal motor 
considered so far; we shall consider some of these differences and attempt to decide 
how they will cause the simple conclusions to be modified. 


5.1. DETERMINATION OF L* 


Consider the V2 motor as an example. The following estimated data will be 
used in the evaluation of equation (14): — 
Sonic velocity: a=3,000 ft./sec. 
Droplet size: r,=0-005 cm. (d,=100 1). 
Vaporisation data: (k/cp,) log. (1 + B)=2 x 10-* cm.?/sec. 
yt1 
Gas dynamic data: [2/(y+1)+(G/pa)? (y—-1)/(v+ 9-91. 
Injection parameter: x,=0:2. 


Hence L*=37°5 &* ft. 
But, taking S=0-5, equation (20) gives 
&k*&=0-14, 
Hence L*=5:25 ft. 
= 63 in. 


The actual value of L* of the V2 motor was 113 in., nearly twice that calculated 
as the minimum necessary. Of course the droplet diameter used in_ thes 
calculations (100) was a guess; if the diameter were doubled, the required L* 
would appear as 252 in. In the absence of data on the droplet size in the V2, ¥ 
can merely conclude that the Model I theory gives results of the right order af 
magnitude. 
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52. IMPORTANCE OF CHEMICAL REACTION RATE 
In the previous section it was assumed that chemical reaction-rate limitations 
were absent. It is necessary to check this. 


To be on the safe side, we demand that E/x, should be less than 0-33, i.e. that 
E should be less than 0-066. Taking the same values as before, and p=0-1 Ib./ft.°, 
it follows that we must have 


(x 10-3)? 
=121 Ib./ft. sec. 


The question is: Can liquid oxygen and ethyl alcohol react together at a 
volumetric rate as high as 121 lb./ft.* sec. at 2,700°K and 220 lb./in.? abs.? The 
data to answer it directly are not available, but the following indirect argument is 
relevant. 


Kerosine and air can react together at 2,000°K and 14-7 lb./in.? abs. at a rate 
of about 75 lb. mixture/ft.* sec. (Ref. 5). Since the reaction rate is proportional to 
pressure to the power 1-8, at 220 Ib./in.? abs., a rate of 10* lb./ft.° sec. would be 
expected. An increase of temperature to 2,700°K would easily increase the reaction 
rate a further fifty-fold. If kerosine and ethyl alcohol are supposed to have the same 
reaction rate, it is clear that the chemical loading parameter E in the V2 motor was 
between one-thousandth and one ten-thousandth of the value which would produce 
extinction. The neglect of chemical limitations in Section 5.1 is therefore fully 
justified. 


5.3. INFLUENCE OF INJECTION VELOCITY 


It has been seen that an increase of x, always increases ¢*, and so the minimum 
combustion chamber volume. It appears therefore that the injection velocity v, 
should be low. In practice this means that the pressure drop across the injector 
orifices should be low, but low injection pressure drop is known to make the rocket 
liable to develop flow instability (“chugging”); moreover fineness of atomisation 
may deteriorate if the impingement velocity of the jets is too low. It is clear that 
2, cannot be made indefinitely small. However, v, is not the only quantity in xo, 
which is defined as pv,/G; it is therefore possible to lower x, by increasing G. 


G is the mass flow rate per unit area in the combustion chamber. It can clearly 
be increased, for given mass flow and chamber volume, by making the rocket long 
and thin, instead of short and fat. Then droplets do not penetrate the flame 
prematurely. Here, too, a limit is set by the difficulty of cooling the walls of such 
acombustion chamber. G also enters the expression for L* to a minor extent; too 
large a value brings, in addition, a thrust penalty. 


Also x, has been seen to enter the chemical criterion which rockets must 
satisfy: a reduction of x,, such as has been advocated, brings extinction nearer. 
\Ithough the V2 example suggested that this is a negligible danger, it must not 
be forgotten that in real rockets the gas in the neighbourhood of the injector, 


February 1959 23 


lated 
hese 
Lt 
we 


D. B. SPALDING 


where extinction would originate, may have a composition and temperature, as q 
result of differing vaporisation rates of fuel and oxidant, which render it much 
less reactive than the stoichiometric mixture. If, in addition, the fuel/oxidan, 
ratio delivered by the injector is not uniform, it is possible that extinction may 
occur locally. This might give rise to unsteady combustion (“screeching”). This 
remark is to be treated as speculative. 


5.4. INFLUENCE OF REYNOLDS NUMBER 


The calculations which have been presented have assumed that the Reynolds 
number of the relative motion of the droplet is sufficiently small for Stokes’s Law 
to be obeyed. This assumption represents an underestimate of both the drag and 
the vaporisation rates. It is shown in Ref. 1 that, when this effect is accounted 
for, the L* value may be halved in a typical case. The foregoing calculations are 
therefore conservative. 


5.5. INFLUENCE OF TEMPERATURE TRANSIENTS IN THE DROPLET 


For simplicity, it hag been assumed that the temperature at which the droplet 
is injected is sufficiently near to the boiling point of the liquid for the transfer 
number to be uniform fhroughout. In practice, however, droplets are usually 
injected at a lower temperature; their vaporisation rate is therefore somewhat 
reduced at first. 


This effect has not yet been investigated in general terms, but it is clear that 
the minimum volume of combustion chamber will be increased over that calcula- 
ted. Priem“ has demonstrated this in particular cases. Further study is needed 
before it can be decided whether the effect is of sufficient magnitude for the motor 
designer to pay attention to it. 


5.6. INFLUENCE OF NON-UNIFORM PROPELLANT SPRAYS 


The assumption that all droplets have the same size and vaporise at the same 
rate on injection is an extreme simplification. In practice the fuel and oxidant 
have different properties, and moreover a wide range of droplet sizes is produced. 
How must our conclusions be modified? 


This matter is currently being investigated. However, some preliminary 
remarks can be made. 


(a) Since the assumption S=0 (zero drag) does not give very poor results, 
and since it is mainly through the drag forces that droplets of one kind influence 
droplets of another, as a first approximation the droplets may be regarded as 
independent of each other. The overall performance of the propellant sprays 
may therefore be calculated by superposing the contributions of the different 
groups of droplets. This means that sufficient chamber volume must be provided 
to burn the largest droplets injected. 


(b) Closer study reveals that the effect of drag renders this conclusion ovét- 
optimistic. Consider the extreme case in which one (uniform) group of droplets 
vaporises and burns immediately on injection, while the second (uniform) group of 
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droplets controls the length of chamber. Instead of the gas velocity being zero at 
the injector face it is now finite. The second group of droplets therefore travels 
faster than the Model Ia theory would suggest and requires a greater length of 
chamber to complete their combustion. A simple calculation leads to the follow- 
ing formula for £*, replacing equation (20):— 


Xo +435 (1 — 0°42) 
2+S 


(22) 


where « is the mass proportion of propellant in the second group, and the refer- 
ence droplet properties are those of the second group. 


(c) Calculation of L* is therefore best made by evaluating ¢* for the largest 
droplets from equation (22), putting «=0. This is conservative. 


5.7. INFLUENCE OF THE FORMATION OF MONO-PROPELLANT DROPLETS 


Many injectors cause jets of fuel and oxidant to impinge together; the subse- 
quent spray therefore probably contains some droplets containing both fuel and 
oxidant. These droplets may be expected to burn somewhat like droplets of 
liquid mono-propellant. It is important to inquire whether the burning rate is 
increased, and if sv by how much. 


This matter has been studied in Ref. 7, where it is shown that the rate is 
indeed increased, the amount depending on the laminar flame speed of the mono- 
propellant and on the droplet size. The effect is greatest for large flame speeds 
and large droplets. Reference 7 contains formulae from which the burning rate 
can be evaluated in particular cases. Application of these formulae to combustion 
in rockets is currently being made. 


5.8. INFLUENCE OF OPERATION ABOVE THE CRITICAL POINT 


This section is concerned with one of the most serious objections which can 
be levelled at the theory presented in this paper: many rocket motors now operate 
at so high a pressure that the distinction between liquid and vapour is absent. i.e. 
at pressures above the critical point of one or both propellants. At such pressures 
there is, for example, no such thing as a latent heat of vaporisation, so the transfer 
numbers cannot be evaluated. Are any of the foregoing conclusions valid for 
these cases? 


This matter has not been subjected to close examination, but the following 
temarks can be made. 


(a) If only one propellant is above its critical point, the present theory is 


applicable to the other one. The same sort of modification is necessary as was 
mentioned in connection with non-uniform sprays (Section 5.6). 


(b) The burning and mixing of jets of fluid above their critical point is much 
more akin to combustion in diffusion flames than that of droplets. Some studies 
of these have been made, and the main principles are clear. However, the 
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theoretical problem is much more difficult than that studied here because partial 
differential equations have to be solved and because the temperature dependence 
of viscosity, density and conductivity must be allowed for. 


(c) In the absence of detailed studies, we can say that, for given x,, the 
length of chamber needed for complete combustion will be proportional to 
d’v,pi/u, where d is the thickness of propellant sheet flowing away from the jet 
impingement zone, v, is the injection velocity, p,.is the density of the injected 
propellant, and « is the gas viscosity. Thus d takes the place of r,; fine 
subdivision of the fuel is still essential. 


(d) Below, but near, the critical pressure, the droplet burning theory is 
anyway only approximate, for the latent heat of vaporisation then ceases to he of 
dominant importance. 


Summing up, although it seems likely that the theory of the present paper is 
qualitatively correct also for rockets operating above the critical point, this matter 
requires urgent study. 


6. Conclusions 
Some of the conclusions which have emerged from this study and which bear 
on the motor designer’s task are as follows: — 


(a) The atomisation of the propellants must be fine. L* is proportional to 
the square of the diameter of the large droplets. 


(b) The injection velocity should be as small as is consistent with absence 
of instability and good atomisation. 


(c) Pre-heating the propellants and forming droplets containing fuel and 
oxidant simultaneously are advantageous. 


(d) In bi-propellant motors, such as the V2, and still more those operating 
at higher pressures, chemical reaction rate is not limiting, except possibly 
in local regions close to the injector. 


(e) Long thin chambers tend to have smaller L* than short fat ones, but are 
more difficult to cool. 


(f) The gas velocity in the combustion chamber increases steadily with 
distance; chamber heat transfer should therefore be calculated from 
stagnation-point data rather than pipe-flow formulae. 


(g) Near and above the critical point, it is more useful to use diffusion- 
flame models than droplet-burning ones. Only qualitative information 
is available as yet. 


Some of these conclusions may, in retrospect, appear obvious. However the 
provision of a mathematical theory for them may be thought to have strengthened 
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their acceptability and given them quantitative meaning. This theory only repre- 
sents a first step; if backed up by careful observations of rockets and the collection 
of the appropriate thermodynamic and transport data for the propellants, more 
useful results can be expected. 
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The Ground Effect on the Jet Flap 
in Two Dimensions" 


D. J. HUGGETT 


(Department of Aeronautical Engineering, University of Southampton) 


SumMaRY: The effect of the ground on a jet-flap aerofoil in two-dimensional 
flow is investigated, with particular reference to the conditions under which the 
jet flow hits the ground. Experimental results are included for the change of 
lift and pitching moment for 58-1° and 31-4° jet-deflection angles, with a 
limited investigation into the downwash changes behind the model with a 
58-1° jet-deflection angle. It is concluded that there is a definite limit to the 
jet coefficient, which may be used at any ground position, if the problems of a 
sudden change of lift and pitching moment are to be avoided. There remains 
the unresolved difficulty of the downwash variation at the tail arising from 
the presence of the ground. 

A provisional assessment of the performance of a jet-flap aircraft near the 
ground is made, and shows that for a 58-1° jet-deflection angle there is a 
definite limit to the minimum flying speed for any given distance of the 
wing from the ground. 

The fact that there appears to be a maximum pressure lift coefficient, 
independent of the jet-deflection angle, which may be obtained at any ground 
position, is investigated. An attempt is made to provide a theoretical 
explanation of this on the basis of the critical condition of blockage underneath 
the aerofoil. 


1. Introduction 


Considerable work, both theoretical and experimental, has been done over the 
past few years on the jet flap. In assessing the scope of this work it was felt that 
little had been done on one particular aspect—the ground effect. The main advan: 
tages of the jet flap must be the reduction, for an aircraft of given weight, of one 
or more of the following quantities: wing area, landing run, take-off run, engine 
size. The reduction of any one or more of these factors requires that the lift 
coefficient at take-off shall be accordingly higher. In this paper an investigation 
has been carried out into the effect of the ground on a jet-flap aerofoil at high 
lift coefficients. 


The basic models used in these investigations were constructed at the National 
Gas Turbine Establishment. Each of the models has a chord of 8 in. and a span 


*This work forms part of the author’s Ph.D. thesis, to be submitted in May 1959. 
Received July 1958, 
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of 12 in. Two models were available, the first with a 58-1° jet-deflection angle 
and the second with a 31-4° jet-deflection angle. The model under test was mounted 
vertically between two horizontal endplates and tested in a 5 ft.x5 ft. return 
circuit type wind tunnel, with a working section length of 7 ft. 6 in. and a speed 
range 0-100 ft./sec. The size of the endplates was chosen in accordance with the 
suggestions given in Ref. 1. Because of the nature of the model mounting (i.e. 
between finite endplates) the experimental set-up does not represent true two- 
dimensional conditions. A calculation based on Fig. 2(a) of Ref. 2 showed that 
the effective aspect ratio of the model was 9-6. Since, however, the lift and pitch- 
ing moment were obtained from chordwise pressure tappings at the mid-span 
station, the results should be closely analogous to two-dimensional results. 


The ground effect was simulated by introducing a movable ground board which 
completely spanned the endplates. The jet coefficient was varied by changing the 
pressure of the compressed air supply to the model. There was a limit to this 
variation, since for structural reasons the maximum internal pressure used was 
15 lb./in.? gauge. Further alteration of the jet coefficient was achieved by reducing 
the main stream speed. 


NOTATION 
c chord of jet-flap aerofoil 


C; jet coefficient, =J/(qS) 
C,, total lift coefficient, =C,,+C, sin 6 
Cy» pressure lift coefficient 
Cy total pitching-moment coefficient about the mid-chord point 
C, pressure coefficient 
d_ vertical distance of aerofoil chord from ground 
J total jet reaction or momentum flux at the nozzle 
I distance of the centre of lift behind the mid-chord point 


M_ net downstream mass flow of the main stream between the aerofoil 
and the ground 


q dynamic head of the undisturbed main stream, =4pU? 
S plan area of jet-flap aerofoil 
U_ undisturbed main stream velocity 
Up,Uag net velocity at P, AB 
w wing loading of hypothetical jet-flap aircraft 
x,z Cartesian co-ordinates defined by Fig. 20 


T’ total vorticity on the aerofoil 
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Ficure 1. The ground effect on the lift coefficient for a 31-4° flap. 
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JET COEFFICIENT om 
Ficure 2. The ground effect on the lift coefficient for a 58°1° flap. 


2 The Ground Effect on the Lift Coefficient 


Earlier experimental work": “ has indicated that there is a loss of lift due to 
the presence of the ground. In the present investigation the lift results (Figs. 1 and 
2) show that the effect of the ground on the lift coefficient is negligible below a 
certain critical jet coefficient. Above this value, however, there is a marked loss of 
lift due to the ground effect. For example, with the 58-1° model at d/c=0-75, the 
critical jet coefficient is 1:25, and at a C, of 3-0 the pressure lift is only 56 per cent 
of the pressure lift at d/c=2-5 at the same C,. 


The critical jet coefficient, at which this loss of lift occurs, varies with the 
position of the ground and jet-deflection angle, and seems to be closely related to 
the jet coefficient at which the jet flow first hits the ground. When the jet flow hits 
the ground, a part of the jet flow is turned upstream and the pressure distribution 
on the ground is altered, the place at which the jet hits the ground being clearly 
defined by a positive pressure peak. By flow visualisation and pressure plotting 
the ground, the critical value of the jet coefficient may be found. Comparing this 
Value with that obtained from Figs. 1 and 2 indicates the close relationship between 
the conditions at which the loss of lift first occurs and the conditions at which the 
jet flow first hits the ground. 


February 1959 31 


12 

Ie 
| | | | 
d 


D. J. HUGGETT 


It is interesting to note that the maximum pressure lift coefficient (C,,) obtain. 
able at any ground position seems to be the same for each deflection angle (see the 
experimental points in Fig. 21). 


3. Ground Effect on the Aerofoil Pressure Distribution 


The pressure distribution diagrams for the aerofoil (Figs. 3-6) have been 
included to show exactly where on the aerofoil surface the loss of pressure lift 
occurs for the 58-1° model. Figure 3 shows a pressure distribution diagram when 
no loss of lift has occurred. An increase in C; now causes a marked change in 
shape, particularly over the rear 50 per cent of the lower surface. Here the positive 
pressure falls and, at higher values of C;, becomes a suction pressure; this loss of 
lift near the trailing edge will give rise to an increment of pitching moment. It can 
be seen clearly in Fig. 6 exactly where this loss of lift occurs and, as well as the 
severe suction peak near the trailing edge, the suction pressure over the whole of the 
upper surface has been reduced. At C;=2-62, so severe is the loss of pressure lift 
that over the rear 40 per cent of the chord the net lift is now zero. For a 31-4° 
flap the loss of lift occurs in a similar manner, but since the loss is not so severe the 
effect on the pressure distribution is not nearly so marked (see Reef. 5). 


4. Pressure Distribution on the Ground 


A typical example of the pressure distribution on the ground is shown in Fig. 7; 
further examples may be seen in Refs. 5 and 6. As mentioned in Section 2, the 
conditions at which the jet flow first hits the ground are clearly seen by reference to 
the pressure distribution diagrams. For, where the jet flow strikes the ground, 
there is a positive pressure peak, the maximum pressure coefficient of which is 
greater than one. When the jet flow first hits the ground the pressure peak is small, 
but, as the C, increases, the pressure peak becomes larger, with values of C, much 
greater than one. It is interesting to note that the point at which the jet flow hits 
the ground moves upstream as C, increases. This, at first, might seem obvious, 
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Ficure 3. Pressure distribution around the aerofoil. C,=0-505. 
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FiGureE 4. Pressure distribution around the aerofoil, C, =1-04. 
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Ficure 5. Pressure distribution around the aerofoil. C, = 1-525. 
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Ficure 6. Pressure distribution around the aerofoil. C, =2°62. 
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FIGURE 7. Pressure distribution on the ground. (Note: The zero for the horizontal scale is the 
point on the ground directly below the model trailing edge.) 


since it would be expected that, as C, increased, the jet flow would tend to approach 
the shape it would have if it were undeflected from its original angle of ejection 
This would certainly be expected over the first few inches of the jet flow. It is 
therefore important to note that, not only does the jet flow approach its “undeflected fF Ficu 
shape” but, with increasing C;, it is bent even farther upstream. 
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FicuRE 8. Flow pattern downstream of the mid-chord point (@=58:1°, d/c=0°5, C,=0°505). 
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Figure 9. Flow pattern downstream of the mid-chord point (@=58°1°, d/c=0°5, C,= 1-04). 
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FiGuRE 10. Flow pattern downstream of the mid-chord point (6=58-1°, d/c=0°5, C,=2°62), 


WIND DIRECTION 


TRAILI 


FiGure 11. Flow pattern upstream of the trailing edge (6=58-1°, d/c=0°5, C, =2°62). 
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5. Ground Effect on the Flow Pattern 


Flow visualisation was obtained by painting a black plate with a mixture of 
titanium oxide, oleic acid and paraffin. It was difficult to find the best proportions 
for these constituents, owing to the extremes of velocity occurring over the plate. 
From the photographs it will be seen that, despite careful mixing and painting, the 
solution has, in places, remained almost unaffected, while in other places it has 
tended to “build up”. Figs. 8-10 show the flow pattern downstream of the 
mid-chord point; Fig. 11 shows the flow pattern upstream of the trailing edge. 


The development of the flow mentioned in the previous paragraphs may be 
traced in the photographs (Figs. 8-10). With low values of C,, the jet flow is turned 
to flow downstream parallel to, and not touching, the ground. This is shown in 
Fig. 8, where it will be seen that the jet flow very nearly touches the ground. In- 
creasing C, causes the jet flow first to touch the ground and then, with further 
increase of Cj, causes the jet flow to hit the ground. At a certain value of C,; the 
jet flow will follow its “ undeflected shape ” (see Section 4)—this condition is shown 
in Fig. 9. Much of the jet flow which is turned upstream is entrained back into 
the jet along its lower boundary when the jet flow hits the ground. This forms an 
eddy which is marked by the build-up of paste just above the ground on the under- 
side of the jet flow. Further increase in C, (Fig. 10) causes the jet flow to be bent 
farther upstream, beyond its “ undeflected shape ”, and a greater proportion of it to 
be turned upstream. The eddy on the underside of the jet flow moves, not only 
upstream, but also away from the ground. Fig. 11 shows the flow pattern upstream 
of the trailing edge for the same C, (2°62) as Fig. 10. A part of the jet flow turned 
upstream at the ground has to pass over the leading edge, completely blocking the 
main stream and causing all of it to pass over the aerofoil. 


6. Ground Effect on the Pitching Moment Coefficient 


Pitching moments'*: ’) have been evaluated about the mid-chord point, the point 
which approximates to the centre of pressure lift with no ground. Figs. 12 and 13 
show the pitching moment for each deflection angle, plotted against the ground 
position for given values of C,. It will be seen that the pitching moment is 
negative throughout the range tested; this is due to the very large nose-down 
pitching moment from the jet reaction. The effect of the ground on the pitching 
moment is not very great for the 31:4° model, but there is a marked change with 
the 58-1° model (Fig. 12). With the larger deflection angle there is, as mentioned 
in Section 3, a considerable increment of nose-up pitching moment as the ground 
distance decreases, due to the loss of pressure lift near the trailing edge. 


For any given C, there exists a critical value of d/c, above which the pitching 
moment is unaffected but below which the pitching moment decreases rapidly. This 
titical ground position for a given C, is analogous with the critical jet coefficient 
for a given ground position (see Section 2). 


In view of this abrupt change in the pitching moment for the 58-1° flap, it is 
perhaps a little surprising to note that the centre of total lift varies very little with 
ground position at moderate values of C, (C; < 2-0) (see Fig. 14). 
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Ficure 12. Pitching-moment coefficient for the 58-1° flap. 
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FicureE 13. Pitching-moment coefficient for the 31:4° flap. 
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FiGurE 16. Downwash behind the 58-1° flap. 


%7. Ground Effect on the Downwash 


The flow pictures obtained led to the consideration of the possible effect of the 
presence of the ground on the downwash at a typical tailplane position. Hence, the 
downwash angle «, measured from the chord line, has been determined by the 
insertion of a yawmeter. Values of the downwash angle were measured along a 
vertical line three chords behind the model mid-chord, ranging from one chord 
above the extended chord line to half a chord below. The readings taken at two 
points on this line, for the 58-1° model, are shown in Figs. 15 and 16 and it will be 
seen that the downwash angle at these points varies considerably as the ground 
approaches the model. For instance, at a point on the extended chord line three 
chords behind the model mid-chord, the change in the downwash angle caused by 
bringing the ground in from 2:5c to 1:0c is 6° at C;=1-805 and 10° at C,=3-15. 
It would appear that nowhere within the range tested is there any substantial 
reduction in these changes. 


8. General Comments on the Ground Effect 


This work has so far been confined to two-dimensional (2-D) conditions. It 
would seem reasonable, therefore, as a first approximation to suppose that the 
ground effect in three-dimensional (3-D) conditions has the same overall effect on 
the total lift coefficient as the 2-D ground effect. There are two factors which may 
significantly alter the ground effect in 3-D conditions. The first is that the jet flow 
will penetrate the main stream more easily, since the main stream flowing under 
the wing may now have a spanwise velocity component. Secondly, when the jet 
flow hits the ground, it may spread spanwise as well as upstream and downstream. 
Of these two effects it would appear that the first would make the ground effect 
more severe in 3-D conditions, whereas the second would tend to reduce the ground 
effect in 3-D conditions. It therefore seems reasonable to assume that the effect of 
the ground in 3-D conditions will be similar to that for 2-D conditions, especially for 
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Ficure 17. Minimum flying speed for an aircraft with a 31:4° flap. 


high aspect ratios. In view of this is seems fair to apply the 2-D results to a finite 
wing, to give an indication of how the ground will affect the performance of a 
hypothetical jet-flap aircraft at low speeds near the ground. 


This performance is shown for an aircraft of given wing loading (w=71-5 
lb./ft.2) in Figs. 17 and 18, where the minimum flying speed is shown plotted 
against wing position for varying thrust/weight ratios, for both 31:4° and 58-1° 
jet deflection angles. It can be seen that, for low thrust/weight ratios, the effect 
of the ground is only significant when the wing is very close to the ground. If. 
however, a 50 per cent thrust/weight ratio is chosen, the ground effect is apparent 
(for the 58-1° jet flap) when the wing is two chords from the ground. From Figs. 17 
and 18 it would not appear to be a practical proposition to use a 31:4° jet flap for 
landing or take-off, since, to achieve the same minimum flying speed as its 58-1° 
counterpart, between two and three times as much power is required. 


For an aircraft with a 58-1° jet flap, Fig. 18 shows that, for a given wing 
position, the rate of decrease of the minimum flying speed tends to zero with increas- 
ing thrust / weight ratio, i.e. a point of extrema for minimum flying speed is obtained 
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FiGuRE 18. Minimum flying speed for an aircraft with a 58-1° flap. 


at a given thrust/weight ratio. This value of the thrust/weight ratio is that which 
gives rise to the critical jet coefficient. It would seem reasonable, therefore, to limit 
the maximum jet coefficient used in practice to that of the critical jet coefficient 
corresponding to the given wing position. This limitation avoids any problems 
that might arise from a rapid loss of lift or change in pitching moment associated 
with the performance of a jet-flap aircraft flying at a jet coefficient greater than the 
critical value. Reference to Fig. 18 will show that, for reasonable wing positions, 
the thrust/weight ratio required (under the given limitation) is small; for instance, 
if the wing position is one chord from the ground, the thrust/weight ratio required 
is 334 per cent. 


There is still the problem of the change in downwash due to the presence of the 
ground, which is not eliminated by this limitation on C,. It would seem that this 
aspect of the suitability of the jet-flap aircraft should be carefully considered, since 
it is possible that this change in the downwash angle due to the ground effect will 
occur in an extensive region behind the wing. It also seems reasonable to assume 
that a canard configuration should be evaluated, in the light of possible upwash 
changes due to the ground effect. 


9. An Approach to the Theoretical Study of the Ground Effect 


Attempts have been made in the past'*:*) to provide a theoretical explanation 
for the ground effect on a jet flap, using a single lifting line as a model. The effect 
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FicuRE 19. The mixing of the jet and the main stream. 


of the ground has been determined as the change of velocity at the lifting line due 
to the induced velocity from its image. 


While, for a first indication of the ground effect, the model would appear to be 
a reasonable approximation, it is clear that there is a marked difference between the 
theoretical values obtained and the experimental values (see, for example, Fig. 1 of 
Ref. 2). This difference may be due to the fact that what may prove to be a 
significant feature has been overlooked. This feature, it is suggested, is the block- 
age of the main stream flow on the underside of the aerofoil, as defined in 
Section 10. If the conditions under which this blockage first occurs are taken to be 
the limiting factor in the increase of pressure lift, it is then necessary to seek an 
expression for the mass flow between the aerofoil and the ground, and to determine 
the limitation imposed on the pressure lift by this condition of blockage. 


10. Blockage 

To determine more exactly what is meant here by the term “ blockage”, a jet 
ejected from a slot AB in a wall XX’ (Fig. 19) is considered. Mixing occurs between 
the jet and the main stream along both boundaries of the jet and gives rise to the 
mixing regions CAD and EBF. If PP’ is a streamline of the jet, then AQ may be 
drawn such that 


mass flow across AP=mass flow across QP’ 
and hence there is zero net mass flow across AQ. 


“ Blockage” between XX’ and YY’ is here defined as that condition under 
which Q becomes a point on YY’. AQ may then be called the blockage line and, 
although the net mass flow across AQ is zero, there is an interchange of main 
stream and jet across this line. It is also true, of course, that the net mass flow 
across a line drawn from any point on XA to any point on YY’ will be zero under 
this condiu..a of blockage. 
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___GROUND 


Figure 20. The theoretical model. 


11. Maximum Pressure Lift Coefficient in the Presence of the Ground 


Consider, in two-dimensional flow, a lifting aerofoil which is represented by a 
uniform distribution of vorticity along the chord. (This is not unreasonable for a 
jet-flap wing, especially at low values of C; —see Fig. 8(b) of Ref. 3; it is, inciden- 
tally, a more accurate representation than a lifting line, particularly at low values of 
d/c.) Let the vortex density be I'/c per unit length along a chord c, distance d 
above the ground, the effect of the ground being equivalent to an image distribution 
A’B’, distance d below the ground (see Fig. 20). 


The net downstream velocity, U», at P due to vortex distributions AB and A’B’ 
in a uniform stream U is given by 


> 
Up=U- (1) 
TC 


2z 2 (2d-—2z) 


Consider the net mass flow, M. between the aerofoil and the ground, where 
d 
M = { Up dz. 
0 
Under the condition of blockage this net mass flow is zero, and thus 
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Ficure 21. Maximum pressure lift coefficient for a jet flap near the ground. 


The lift per unit span may be written as pU,,l', where U4, is taken as the velocity 
at the mid-chord point, which, from (1), is 


Cc 


Hence the lift per unit span becomes 


p (v tan 4d Ci pU?c, 
where C,, is the pressure lift coefficient. Rearrangement gives 


2 


From equations (2) and (3), the pressure lift coefficient at the condition of blockage 
may be found in terms of d/c. In Fig. 21 this curve is plotted with the experimental 
values taken from Figs. 1 and 2, and there is very good agreement. 
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12. Conclusions 


Although the effect of the ground is far less on a 31-4° jet flap than on a 
58:1° jet flap, there appears to be little advantage in using a 31-4° jet flap since, 
to achieve the same minimum flying speed as its 58-1° counterpart, between two 
and three times as much power is required. 


It seems reasonable, for a 58-1° jet flap, to limit the maximum jet coefficient 
used in practice to the critical jet coefficient. This avoids the problems associated 
with the loss of lift and change of pitching moment as the wing approaches the 
ground, but does not solve the difficulty of the change in downwash behind the wing. 


At any given wing position there appears to be a maximum pressure lift. This 
value seems to be independent of the jet-deflection angle and may be obtained 
theoretically as dependent on the condition of blockage between the wing and 
the ground. 
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The Vibrations of a Thin-Walled Stiffened 
Cylinder in an Acoustic Field 


J. H. FOXWELL, B.Sc.(Eng.) 
and 
R. E. FRANKLIN, B.Sc.(Eng.), Ph.D. 


(Department of Aeronautical Engineering, University of Southampton) 


SUMMARY: When a vibrating structure encloses a volume of fluid, the acoustic 
effects within this volume modify considerably the response characteristics of 
the structure, provided that the cylinder is vibrating in radial modes only. 
Measurements made of the displacement caused by a particular sound wave 
are of the same order as the values predicted. The calculation of the response 
of the cylinder to an acoustic wave also yields the sound field inside the cylinder 
and, again, the results are in general agreement with practical experience. 


1. Introduction 

During the past few years the problem of vibrations induced in aircraft 
structures by noise has steadily grown in importance. There are two main reasons 
for this; in military aircraft the vibrations set up have been serious enough to cause 
secondary structural failures and failures of electronic equipment mounted within 
the fuselage, and in civilian aircraft the vibrations allow sound to be transmitted to 
the passenger cabin to a degree which affects the comfort of the passengers. It is 
becoming clear that the problem is one which needs to be considered in the design 
stage of an aircraft. The information required is, in the first case, a description of 
the vibrations of the structure suitable for fatigue estimations and, in the second. 
an estimate of the level to which the noise in the fuselage will build up. The basic 
problem is the same in both cases, the only difference being one of emphasis. 


The overall problem has been discussed by Powell who suggested that the 
problem be divided up in the following way : — 


(a) Description of the pressure field causing the vibrations. 

(b) Determination of the modes of vibration of the structure. 

(c) Calculation of the effects of acoustic damping and acoustic reactance. 
(d) Determination of the structural damping. 


The present paper neglects the structural damping and, further, does not 
concern itself with the description of the incident pressure field which, with a noise 
field, is a formidable problem in itself. 


Received March 1958. 
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The main concern of this paper, then, is the calculation of the modes of 
vibration of a structure taking into account the surrounding medium. Hitherto, 
most investigations into the vibrations of structures have ignored the presence of 
the surrounding medium. It has been assumed that, since the disturbances are 
acoustic and therefore small, they will not have an appreciable effect on the motion 
of the structure. However, this assumption is not necessarily true and should be 
questioned, particularly when the vibrations are in fact induced by acoustic 
pressures or when the structure under consideration is very lightly damped and 
probably sensitive to energy losses involved with the radiation of sound. This 
problem of acoustic effects in structural vibrations has been discussed in an 
elementary way in Ref. 2, and it was demonstrated that under certain conditions 
they become most important. Such conditions are met when the structure encloses 
a volume of fluid, a good example of this being an aircraft fuselage. 


The formulation of the problem dealt with in this paper has been influenced 
by the questions of whether the overall problem may, in fact, be divided as just 
suggested and whether reflection and scattering effects are important. As a result 
a simple and relatively idealised problem has been chosen which can be solved 
exactly: it is the two-dimensional problem of an infinitely long stiffened cylinder 
in a field of plane simple harmonic sound waves, the wave fronts of which lie ina 
plane which is parallel to the plane formed by the cylinder axis and a cylinder 
generator*. 


In this two-dimensional problem the cylinder is subject to radial and tangential 
movements only, there being no modal deformation in the axial direction. Here 
the displacements may be represented by a Fourier series, each term corresponding 
to a mode of deformation and the coefficients to the generalised Lagrangian 
co-ordinates. The cylinder is subject to an incident acoustic field and scattered 
and internal fields, the effect of which is to produce generalised forces Q, for each 
mode n. The total kinetic energies (J,) and strain energies (S,) for the shell and 
stiffening members may be written in terms of the generalised co-ordinates. These 
have been given by Miller, who has extended the work on uniform shells‘ ” to 
stiffened shells. The Lagrangian equations may now be used :— 


dt 


=Q,. 


oT: OS: 
) 


For each mode, four equations result, the solution of which gives the values of the 
generalised co-ordinates and hence specifies the cylinder displacement under the 
action of the acoustic fields. 


*It is hoped that more complicated problems will be dealt with in future papers. 
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f NOTATION 
), a radius of cylinder 


r,@ co-ordinates (see Fig. 1) 
n z_ local radial distance 


w radial displacement from equilibrium 


d v tangential displacement from equilibrium 
is Dns Cn generalised co-ordinates in mode n (see equation (1)) 
. n mode number, =number of full waves round the cylinder 
w frequency 

c speed of sound 
k=o/c 
; p density of air 
An, Bry Cn, Dn Fourier constants 
et ps sound pressure in scattered field 
“ p; sound pressure in internal field 

p. sound pressure in incident field 

Qin, etc. generalised forces in a,-direction 

ial J, (kr), (kr) solutions of Bessel’s equation 
“ H,® (kr) the Hankel function, (kr)—i¥. (kr) 
an R, acoustic damping 
ed X, acoustic reactance 
Xy mechanical reactance 
ese P, amplitude of pressure in plane sound wave 
2h cylinder wall thickness 

ps Shell density 

E Young’s modulus 

Poisson’s ratio 

A,B,C,D_ see equations (38)}(41) 

es A; area of frame 
the I, moment of inertia of frame about skin 


Z distance of frame centre of area from skin centre line 
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@ 
0,008 np+b, SiN np 


SIN np+d COs np 


FIGURE 1. ‘ 
Co-ordinate system. 


2. The Acoustic Fields and Generalised Forces 


2.1. THE ACOUSTIC FIELDS 


Consider a cylinder of radius a, centre O, the origin of co-ordinates r, 9 (Fig. 1). 
The radial and tangential linear displacements of the cylinder surface from its 
equilibrium position are denoted by w, v respectively, with 


The coefficients. which are functions of time only, are the generalised 
co-ordinates whose values are required. 


The cylinder is excited by an incident acoustic field p,(r, 9, t). On incidence 
there is a scattering of sound about the cylinder and a transmission to the interior; 
these two fields are denoted by p,(r.¢, 0), pi(r,¢.0. The scattered and incident 
field pressures combine linearly to give the pressure at a point r.¢(r >a), 


D(r, >, (r, + Ds (1, 0), (2) 


while within the shell there is a field pressure 


Since these are acoustic fluctuations they obey the wave equation 
1 0? 
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The general solution of this equation is, in polar co-ordinates, 


cos no J, (kr) 


exp {iwt}, . (5) 
Y, (kr) 


sin no 
where J, and Y, are the solutions of Bessel’s equation and p is exactly determined 
by the boundary conditions. 


The acoustic fields are now considered in turn. 


(a) The Incident Field 


In a practical case the incident acoustic field will often be a random one, that 
isa noise field, the mathematical representation of which is difficult and requires 
separate consideration. Noise fields may originate from turbulence in a boundary 
layer or in a jet efflux. The nature of these fields is unknown at present. Even if 


it were known, the application to this problem would present considerable 
difficulties. 


For a non-random exciting field the problem is much simpler; for instance, for 
a plane sound wave the incident pressure field is defined (using the co-ordinate 
system of Fig. 1) by 


p=P, exp {ik (x+ct)} =P, exp {ik (r cos ¢+c?)} 


will, (J, (kr) +2 (i"Jn (kr) cosng)exp {iot}. © 


The Scattered Field 


The scattered field is defined, through the boundary conditions, by the geometry 
of the shell and its dynamic deflection. 


The boundary conditions to be satisfied are 


(i) The particle velocity of the scattered plus incident fields at the surface of 
the cylinder (r=a) must equal the velocity of the surface (w). 


(ii) At infinity the scattered wave must have the form of a diverging wave. 
(The Sommerfield radiation condition). 


Condition (ii) is satisfied by the solution to the wave equation (equation (5)) 


(A, cos n¢ +B, sin n¢) H,,© (kr) exp {iwt} (7) 
where H,, (kr) is the Hankel function, =J, (kr) —iY, (kr), and 
Lim H, (kr) exp {iwt} = / {2/(tkr)} exp { (r—q/k-ct)}, 
Which has the form of a diverging wave. 


Condition (i) determines the constants A,, B, in equation (7). 
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(c) The Internal Field 


Again the internal field is defined by the shell and its dynamic deflection. 
The boundary conditions in this case are 


(i) The particle velocity in the internal field at the cylinder surface (r=q) 
must equal the velocity of the surface. 


(ii) At the origin (r=0) the pressure must be finite. 


Condition (ii) leads to the rejection from the solution (equation (5)) of the function 
Y, (kr) which is infinite at r=0. 


The solution is therefore 


=X(C, cos no+D, sin no) J, (kr) exp {iwt}, . 


Mea 


where C, and D,, are determined by condition (i). 


2.2. THE GENERALISED FORCES 


The generalised force Q, is defined in the Lagrangian equations as 


Q . . . (9) 


where W is the work done through the displacement of the loading. The loading 
per unit area arises from the acoustic pressures and is, in the outward radial 
direction, 


The forces acting in the tangential (v) direction are viscous forces, which are very 
small compared with the acoustic normal forces; therefore the generalised forces in 
the c, and d, directions are neglected. 


22 
Now sw= (Sa, cos n¢+5b, .  . ll) 


0 


and the generalised forces per unit length are 


2n 

0 
3 2n 

Q, = =a | {pi—(Po+Ps)} rea inno dd. . (1) 
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On substituting for p,, p, from equations (7) and (8), these become 


2n 


Q., a[ (ka)— | rae 008 n9 do (ka) ] exp {ior} (14) 


= (1) 
2s 
[ (ka) - | {Do} rma Sin n dp (ka) exp {it} . (15) 
ion 
and A,, B., C,, D, are found from the boundary condition that the particle velocity 
and cylinder surface velocity be equal at r=a. The particle velocity (Ww) in a 
sound wave is given by the equation of motion 
(8) _ Op _ 
which, for w=f (r, ¢) exp {iwt}, reduces to 
Application of the boundary condition gives 
(9) ri 
w= a, cos no+ b, sin nd= — 
Jing p 
dial 


1 a 
jwp or {Pot Ds} 


(10) } whence 


very 22 
A,= exp {iwt} (ka) [a nlpC + { or Cos np do | (18) 
0 
2n 
B,=- Sap (ka) [ { sin no dp (19) 
ipc 
(12) 
be 21 


(13) Finally, on substituting these coefficients in equations (14) and (15), the generalised 
forces are found to be as given by equations (22) and (23). 
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rag @) 

Q., inapc — - an +F,, (t) . (22) 
JI, 


2n 


where F, {Po} r=a COS NP { cos n¢ do 
r=a 


(2) 
G, (j= 2 | sin no do sin no do. 
=a 


The first terms in (22) and (23) can be written in a more convenient form by 
making use of the following properties of the Bessel functions : — 


H,©=J,—iY, 
LOY 
where, in this case, x= ka. 
Hence Q., =[RatiX,] dnt+F, (24) 


1 2 2 
where R, = (2) 


1 2pc\ 2 
and X,= J, yale 


R, and X, are the acoustic radiation and reaction terms. Ry, is real and represents 
the damping of the cylinder vibration due to re-radiation of sound waves. 
R,+iX, together represent the acoustic impedance of the system at r=a. 


3. The Lagrange Equations 


3.1. INTRODUCTION 


If the potential and kinetic energies of the cylinder are known in mode 1, 
together with the generalised forces in that mode, then, by means of the Lagrange 
equation, the differential equations for the generalised co-ordinates a,, b,, Crs dr Cal 
be written; the equations can then be solved, giving the disp!acements in mode r. 
These displacements may then be summed (vectorially) over the m modes to give the 
total displacements w and v. 
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Expressions for the potential (strain) and kinetic energies of a uniform shell are 
given by Warburton and Arnold“, who used expressions derived by Love'” for the 
stresses in a Shell element. A full discussion of the derivation of these formulae 
and the validity of the assumptions involved is given by Miller®. Méiller also 
derives more exact expressions and extends the analysis to include the effect of 
stiffening by frames and longerons. For the present calculation the more approxi- 
mate formulae of Ref. 6 are used and are extended to include frame stiffening. 


32. THe POTENTIAL AND STRAIN ENERGIES 


Using the co-ordinate system defined by Fig. 1, the expressions for the energies 
in two dimensions (no axial modes) are (for the uniform shell) 


S,=strain energy per unit length 


0 -h 
where w is positive outwards. 
2h is the cylinder wall thickness (note that the origin of co-ordinates (v, w) is 
the shell middle surface). 
T,, the kinetic energy per unit length is equal to 


22 +h 


0 -h 
where p, is the shell density, and dots denote differentiation with respect to time. 


Substituting for w and wv (equation (1)) in equations (26) and (27) and 
integrating, 


0 En 


In this two-dimensional problem the shell is stiffened by frames only. This is 
done in the interests of relative simplicity and the avoidance of complicated 
‘pressions which result from the inclusion of longeron stiffening. To calculate the 
additional energy due to the frames the integrals (26) and (27) are used, but the 
tintegration now extends over the depth of the frame. Provided that the ratio 
(epth of frame)/(cylinder radius) is small, it is felt that these will give satisfactory 
ipproximations. Since there are no modes in the axial direction, it is assumed that 
the “ frame energies ” are uniformly distributed along the length of the shell. 
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If d is the frame spacing, then from (26) 


S;=frame strain energy per unit length 


2n 


0 depth 


where 5; is the width of the frame (i.e. in the axial direction) and the term Up a? is 
neglected in comparison with v, /a. 


Substituting for w and v and integrating 


Ex At Zz 
2 2 4 Ards 


where | b, z dz, 


depth 


2;= distance from the skin centre line to the centre of area, 


b; dz=second moment about the skin centre line. 
depth 


The kinetic energy is given by equation (27) (with z limits changed). 


T;=frame kinetic energy per unit length 


0 


3.3. THE EQUATIONS OF THE CYLINDER MOTION 


Having found the generalised forces (equations (24) and (25)) and the potential 
and kinetic energies of the shell (equations (28), (31) and (29), (32)), the Lagrangian 
equation 


may be solved. The solution is given by the following equations: — 


Alin + Ban+(Rat+iXa) Gn+CCn =F» (t) . 3 
Aén+De,+Can =0 « 
Ady+Dd,+Cby 
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where A= 2p,ahe (1+ 45). 


The differential equations (34)-(37) form a set of four equations for each value 
of n. Each set is independent and not coupled to any other set; there is no 
coupling of any sort between modes for different values of n. 


These equations can be solved for @,,, bn, C», d, when F, (¢) and G, (t) are known. 
To illustrate the solution, take a simple harmonic exciting field, i.e. 


F,, (t\)=F, exp {it} 
(t)=G,, exp {iwt}. 
Putting a, (t)=a, exp {iwt}, 


F,, 


(Xu — OX) + (42) 
d,= 5 « @& 
where u= —Aw?+ B+ (46) 


Xy may be termed the mechanical reactance of the vibrating shell. In the absence 
of the acoustic reactance (X,) resonance occurs when Xy=0. It will be noticed that 
the coefficients B, C, D are functions of n, the mode number. Thus there is one 
Xy for each mode and the zeros give the resonant frequencies (in vacuo) of the 
shell in each mode. However, for vibration in an acoustic medium, it can be seen 
(equations (42)-(45)) that the total reactance is modified by the acoustic reactance 
X,. Resonance now occurs when 


Xx =0. 
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Examination of X, (equations (24)-(25)) reveals that, at the roots of J,’ (wa/o), Zz, 
becomes infinite; thus, if the root of J,’ occurs at a value of » close to the resonant 
frequency of mode n in vacuo, the response curve changes shape. 


The term iwR, in equations (38)-(41) represents the “ acoustic damping ” of the 
shell, i.e. the energy loss due to the radiation of sound. Structural damping would 
add a further term. 


4. The Response of the Cylinder to a Simple Harmonic Plane Wave 


In the preceding two sections the calculation of the response of the cylinder to 
an arbitrary periodic incident field has been described. To illustrate the character. 
istics of the problem the calculation is carried out in the following sections for an 
incident sound field of plane simple harmonic waves. 


The pressure in a plane wave travelling in the —x (Fig. 1) direction is 
expressed by 


Po=P, exp {ik (x+ct)} 
=P, exp {ik (r cos ¢+ct)}. 
This can be written’ 
(kr) +2 & Py (Un (kr) cos nolexp {iot} . (4) 
n=l 
in which the plane wave is now expressed as an infinite number of cylindrical waves. 
Putting p, into equations (22) and (23), 
(2) 
F, ()=F, exp {iot} =2 
G, ()=0. 


Therefore, from equations (48) and (42), the amplitude of the displacement w 
is given by 


2 c 1 
{(Xy—oX,)? + o?R,?}! {(Xy— oX,)* + w?R,?}4 
| bs [=0. 


The computation of | a, | has been carried out for the mode n=10, using 4 
stiffened cylinder, and details are given in the Appendix. 


The response curves calculated are shown in Figs. 2 and 3. 


5. Discussion of the Response Curves 


Figure 2 shows the mechanical reactance (Xy) and the acoustical reactance 
(wX,) plotted to the same scale for the mode n=10. From equation (44) it will be 
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Xy 


Ficure 2. Mechanical and acoustical reactance for mode n= 10. 
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Ficure 3. Response of cylinder in mode n=10 excited by a plane sound wave. 
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seen that resonance occurs when the frequency of the incident wave is such that 
Xu =0. 


‘Therefore the intersections A, B, C of Xy and wX, in Fig. 2 denote resonant points 
of the system. It will be appreciated that there are an infinite number of resonant 
points, since »X, has a quasi-periodic form and there are thus an infinity of 
intersections. 


Figure 3 shows the response curve for the mode n=10; the resonances A, B,C 
correspond to the intersections in Fig. 2. The broken curve shows the response in 
the absence of acoustical effects. These effects are now considered. 


Acoustic damping arises from the fact that radiation of energy away from the 
vibrating systems represents to the structure a damping of its motion; this damping 
is given by equation (24), 


En i 
and the magnitude of the acoustical damping can be determined for any frequency 


and mode from this expression. Junger” gives curves of R, for n=0(1)10 
plotted against wa/c.* 


Another important damping term is the structural damping arising from the 
hysteresis of the material and from riveted joints. The magnitude of the structural 
damping is not precisely known, but it is probable that the structural damping is 
often comparable in size with the acoustic damping. The effect of damping on the 
response will be to diminish peak values and to increase the effective bandwidth 
about the resonant points. 


The acoustic reactance X, is mainly responsible for the changes in the 
response curve from the in vacuo condition (broken curve). There are three conse- 
quences of the effect. Firstly, that, whereas in the in vacuo case there is only one 
resonant frequency, in the present system there are many resonant frequencies. 
(This is due to coupling between the cylinder vibration and the wave motion within 
the cylinder.) Secondly, there are points on the curve at which the response is, 
theoretically, zero. These arise from standing wave patterns which occur in the 
internal field at certain frequencies corresponding to the roots of the equation 


J,’ (24) =0 (see equation (25)). 


Under certain conditions this gives rise to the third consequence when the frequency 
at which standing waves occur coincides with the in vacuo resonance of the cylinder. 
In this case the normal (in vacuo) resonance disappears. 


Equation (49) gives the amplitude of the motion in the n™ mode; to find the 
total amplitude it is necessary to know the time phase of each modal displacement; 


*While Junger’s paper does not consider the transmitted interior field, his function for Ry, 
is identical with the result of this paper, since there is no energy loss in the internal field. 
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TABLE I 


THE RESPONSE IN THE FIRST 20 MODES TO A PLANE WAVE OF AMPLITUDE P, LB./IN.? AND 
FREQUENCY 577 ¢./s. 


n=19 
The radial displacementw= | a, | exp {i (wt+8,)}=20x x P, in. 


n=0 


| a, | and 8, are tabulated for »=577 c./s. 


|/P, |a, 

n (in. x 10-8) n (in. x 10-3) 

0 5-0 172°27 10 26-0 193°10’ 
1 0:3 166°8’ 11 1:92 27°48’ 
2 1-4 352°44’ 12 1-38 105° 40’ 
3 1:8 362°5’ i3 0:99 126° 19’ 
4 | 2-4 |  164°59’ 14 0:69 185° 

5 | 2:0 | 15 0:66 73°19’ 
6 2°5 |  231°20° 16 0-51 |  152°21’ 
7 | 2: 75°42" 17 | 035 246240" 
8 4-0 104°13" 18 | O18 |  348°42° 
9 | 9-0 |  140°58” 19 | 0-08 86°14’ 


this is readily obtained from equations (42) and (48), which give 5,, in degrees, as 


1 
[ tan (X,, o¥,)+ oR. +90 (n+1)| . (50) 
where x=wa/c and the radial displacement 
oo 
w= > | a, | exp {i (wt + 6,)}. ‘ 


Table I shows the amplitudes and phases of the first 20 modes. A summation of 
these gives a total displacement of 20 x 10-* in. for a plane wave of 160 dB. (datum 
pressure =0-0002 dynes/cm.’) intensity and a frequency of 577 c./s. 


6. The Transmitted Sound Field Within the Shell 


From the expressions already derived it is possible to calculate the distribution 
of sound within the cylinder. Equation (8) gives the value of the interior sound 
field p;, where the coefficients C, and D, are given by equations (20) and (21). As 
would be expected, the interior field is dependent upon the radial motion of the 
cylinder surface (w). The sound pressure is also proportional to J, (kr); thus at 
the centre of the cylinder (r=0) the sound pressure depends principally upon the 
mode n=0, since J, (0) ~0 for n=>1. Therefore, at frequencies remote from the 
resonance of the n=0 mode, the sound pressure level at the cylinder centre may be 
a minimum. 
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FicurE 4. The sound field along ¢=0 for an incident plane wave of amplitude P, and 


frequency 577 c./s. 


Figure 4 shows the sound pressure level along the ¢=0 line for an incident 


plane wave of 577 c./s. 


and thus the sound pressure is high at the centre. 


from the principal mode (n= 10) at this frequency. 


This frequency is close to the resonance of the n=0 mode 


Also plotted is the contribution 


From equations (7) and (18) it is possible to compute the external field pressure 
(p, +p). The expression for p, contains two terms, the second of which represents 
the scattering when the cylinder is rigid and the first the effect of the motion of the 
cylinder surface. Wiener’ shows the scattering from a rigid cylinder. For 


ka> 10 and ¢=0 there is “ pressure doubling” at the cylinder surface. 


Motion 


of the surface causes alleviation of “pressure doubling” and, in the example 
calculated (Fig. 4), there is a pressure increase of 5-8 dB., as opposed to 6 dB. in the 


rigid case. 


7. Conclusion 


In a previous paper the authors suggested that, when a vibrating structure 
encloses a volume of fluid, the acoustic effects within this volume may be such as to 
modify considerably the response characteristics of the structure. The present 
paper has shown this to be true for the two-dimensional problem of a cylinder 


vibrating in radial modes only. 


The principal results of the paper are contained in Fig. 3, which shows how 
the response curve of a given mode of vibration is changed by the effects of the 
surrounding medium, and Table I, which shows that a wave of a given frequency 
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can excite many modes of oscillation. The calculated displacement of 0-020 in. 
for an incident sound wave of frequency 577 c./s. and intensity 160 dB. is of the 
same order as displacements measured in practice. 


Although the main concern of this paper is the calculation of the response of 
the cylinder to an acoustic wave, it has been demonstrated that the same calculation 
serves to give the sound field inside the cylinder. Again, the results, which are 
shown in Fig. 4, are in general agreement with practical experience. 


The method presented here is to be extended to three dimensions, where the 
cylinder vibrates in longitudinal as well as radial and tangential modes. In this 
case it will be possible to calculate the response and internal acoustic fields arising 
from more complex incident sound fields, such as those produced by propellers. It 
should then be possible to arrange the propeller configuration and phasing to 
produce optimum acoustic conditions within an aircraft fuselage. 
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Appendix 


DETAILS OF THE CALCULATION FOR THE PLANE WAVE EXCITATION 
OF THE STIFFENED CYLINDER 


The cylinder constants were as follows :— 
radius =a=60 in. Young’s modulus = E=10' 1b./in.? 
thickness =2h=4-8 x 10~? in. Poisson’s ratio =o=0-3 
density =p,=0-117 1b./in.* 
The frame section chosen was as shown in Fig. 5, whence 
area of frame=A,=0-153 in.? 
moment of inertia of frame about skin=/7,=0-560 in.* 


distance of frame centre of area from skin=z,=1-5 in. 


SDIN. 
|_0-0361N. FIGURE 5. 


| 
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The calculation was carried out for the mode n=10. Using the values given, the 


constants in the expression for Xy (equatioi. (46)) are 


A=2-688 x 10-* Ib.sec.2/in.2 
B=5-393 x 10* Ib./in.2 
C=2:149 x 10° 1b. /in.2 
D=3-170 x 108 Ib./in.2 


Other constants used were 


speed of sound=c=1120 ft./sec. 
density of air=p=2-38 x 10~-* slugs/ft.$ 


The functions J,(ka), Y,(ka) were obtained from Ref. 9. J,’(ka) and Y,/(ka) were 
computed from the relationships 


J,(ka) =} [J,-, (ka) —IJn4,(ka)] 
Y (ka) =4 Yn4,(ka)]. 


Even for one mode, the computation on a desk calculating machine is considerable. 
A complete analysis over many modes would best be done on a digital computer. 


The three-dimensional problem is of more interest and in this case the computa- 


tion would almost certainly be done by a digital computer, since the number of modes 
involved in a given frequency range may be enormous. 
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Thermal Stresses in Rectangular Plates* 


J. S. PRZEMIENIECKI, Ph.D., Grad.R.Ae.S. 


(Bristol Aircraft Limited) 


SumMMARY: The characteristic functions for beam vibration modes are used 
to derive an approximate solution for the calculation of thermal stresses in 
rectangular isotropic flat plates subjected to arbitrary temperature distributions 
in the plane of the plate and constant temperatures through the plate thickness. 
The thermal stresses are obtained in the form of generalised Fourier expansions 
in terms of the characteristic functions, and their derivatives, representing 
normal modes of vibration of a clamped-clamped beam. Since these functions 
have recently been tabulated, the practical application of this new method to 
the thermoelastic stress analysis of plates presents no difficulty. 


1. Introduction 


An approximate variational method for the calculation of thermal stresses in 
thin rectangular plates of constant thickness has been developed by Heldenfels and 
Roberts”’. In this method the Airy stress function is assumed to be given by 
F=f (x) g(y), where f(x) is the stress function corresponding to the thermal stress 
distribution « in an infinitely long plate and g(y) is obtained from the principle of 
minimum complementary energy. Singer et al calculated thermal stresses in 
plates by reducing the problem to that of determining the end effects of the self- 
equilibrating thermal stress o applied at the short sides of the plate and then 
solving the resulting end traction problem by the method of self-equilibrating 
orthogonal polynomials derived by Horvay. This method, however, is only 
applicable to plates with aspect ratios greater than two and with the temperature 
variation restricted necessarily to one direction. Mendelson and Hirschberg used 
a collocation procedure whereby the partial differential equation for the Airy stress 
function is satisfied everywhere in x but at only a finite number of values of y, and 
the problem is thus reduced to the solution of a set of simultaneous fourth-order 
ordinary differential equations for some unknown functions of x. The main dis- 
advantage of the collocation method is the appreciable amount of computing work 
involved in obtaining solutions to the set of simultaneous differential equations 
whose right hand sides depend on the particular temperature distribution 
considered. 


*This paper is based on Chapter 8 of a Ph.D. thesis approved by the University of London‘); 
it was written while the author was at the Imperial College of Science and Technology, 
London. 
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The solution presented in this paper is derived directly from the governing 
differential equation by expressing the Airy stress function as a generalised Fourier 
expansion in terms of the characteristic functions (eigenfunctions) representing 
normal modes of vibration of a uniform clamped-clamped beam™. In so doing 
the problem is reduced to the solution of a set of simultaneous linear equations for 
the unknown generalised Fourier coefficients in the stress function. Since these 
equations are extremely well-conditioned, they are particularly suitable for 
iterative solutions and, in general, one or two iterations are usually sufficient for 
engineering purposes. The overall accuracy depends on the number of coefficients 
evaluated and on the form of the temperature distribution; violent variations in the 
temperature distribution naturally require a greater number of coefficients. The 
solution is approximate in the sense that in any practical problem a finite aggregate 
of terms is used and, consequently, the differential equation for the Airy stress 
function is satisfied only approximately. All boundary conditions, on the other 
hand, are satisfied exactly. 


Since the normal modes of vibration of a clamped-clamped beam are used to 
represent the Laplacian of the temperature distribution V’T, occurring in the 
differential equation, the solution is therefore restricted to temperature distributions 
for which V?7 and the normal slopes of V?T at all four edges would vanish. This 
implies that, for temperature variation in one direction only, changes in temperature 
near the edges must be linear. Such a restriction, however, introduces a negligible 
modification in the actual temperature distribution and the difference between the 
exact stress distribution and one obtained by the present approximate method is 
insignificant. A numerical example is included to illustrate the practical 
application of the method. 


NOTATION 
x,y rectangular co-ordinates, measured along a pair of edges 
f=x/a 
n=y/b 
a,b plate dimensions 
u,v total displacements in x- and y-directions, respectively 
T=T (x, y), temperature distribution 
temperature at »=4. 
Young’s modulus 
a coefficient of thermal expansion 
v Poisson’s ratio 
K=-2E 
F Airy stress function 
direct stresses 


Tz, Shear stress 
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ng V, characteristic functions (eigenfunctions) representing normal modes 
ier of vibration of a clamped-clamped beam 
ng 8B, characteristic value (eigenvalue) 
ng yr constant 
for 
ase Stress function coefficient 
for bin coefficient defined by equation (24) 
Kronecker delta 
the f(x), g(y) functions of x and y, respectively 
‘he m,n, integers 
ate V*=(07/dx* + 0°/0y*), Laplace operator 
ess 
Vt=(0'/0x' +20'/dx*0y* + 0*/dy'), biharmonic operator 
2. Thermal Stresses and Deformations 

to It is assumed that the properties of the plate material do not vary with tempera- 
the ture and that the temperature gradients in the plate are below a critical value 
wd which might produce thermal buckling. The plate is subjected to an arbitrary 
his temperature distribution T (x, y) and there is no variation in temperature through 
ure | the plate thickness. The problem is thus one of generalised plane stress, and it 
ble can be reduced to finding the Airy stress function F which obeys the biharmonic 
differential equation 

is 


and satisfies the boundary conditions of zcro stresses at the edges. 


The co-ordinates x and y refer to a set of rectangular axes along a pair of 
edges but, for convenience in the subsequent analysis, non-dimensional co-ordinates 
are introduced: €=x/a and =y/b, where a and b are the plate dimensions (see 
Fig. 1). The thermal stresses and total strains are obtained from the Airy stress 
function, using the following expressions : — 


= (o,—vo,)+eT . ‘ ‘ ‘ ‘ (5) 
= (c,—vo,)+ aT. ‘ ‘ (6) 
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y,7 é=x/a 
n=y/b 
y=b 
x=a 
FIGURE 1. 


It will now be assumed that the Airy stress function F is represented by a 
generalised Fourier expansion in the form of a double series 


lI 
Me 


m=1 n=1 


F 


where d,,, are the arbitrary coefficients to be determined later and ®,, (€) and ¥, (n) 
(with mand n=-1, 2, 3,.... ) constitute an infinite set of orthogonal functions 
within the region bounded by a rectangle, so that it will be possible to express any 
arbitrary function in this region as a linear combination of the ‘P,, and V, functions. 
This will ensure that the stress function F is complete in the prescribed region and 
that, by taking a sufficient number of terms in the expansion represented by 
equation (7), any required degree of accuracy may be achieved. 


Substituting equation (7) into equations (2), (3) and (4), results in the following 
expressions for the stresses : — 


m=1 n=1 

m=1 n=1 
| oc oo 

ab m=1 n=1 


where primes denote differentiation of the functions ®,, or VY, with respect to theit 
assigned variables. If necessary, the resulting displacements can be calculated 
using one of the standard methods for the determination of displacements from 4 
known stress function. 
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TABLE I 
VALUES OF AND y,. 


(Reproduced from Ref. 6 by permission of the University of Texas) 


4730041 | 0982502 4 | 14137166 | 1-000002 

2 7853205 ‘1000777 5 17:278760 1-000000 

3 | 10-995608 | 0:999967 r>S (2r+1)2/2 | 1-000000 


3, Solution by Characteristic Functions 


Using equations (8), (9) and (10) the boundary conditions of zero stresses at 
the edges become 


®,,=®,/=0 at £=0 and 1 
(11) 


v,=V,’=0 at 7»=0 and 1 


These boundary conditions can be satisfied by taking ®,, (€) and W,,(») to be the 
characteristic functions (eigenfunctions) representing normal modes of vibration of 
a clamped-clamped beam. These functions are derived from the differential 
equations 


@,,”” Bui Pn = 0 
(12) 
of which the solutions are 
(€)=cosh —cos 8,€ — y, (sinh - sin B,€) (13) 
(n)=cosh £,n-cos (sinh 8,9 — sin 
where 8, is the r** positive root of the transcendental equation 
and (15) 


Numerical values of the characteristic functions represented by equations (13) have 
been tabulated in Ref. 6. The values of 8, and y, for r varying from 1 to 5 are 
given in Table I, while for r > 5 it is sufficiently accurate to evaluate 8, and y, from 


3; r>5 . (16) 


February 1959 69 


(1) 
tions 
any 
ions. 
d by 
wing 
(8) 
a 


I. S. PRZEMIENIECKI 


The characteristic functions “,, (¢) or V, (y) have the following properties 


0 
1 
0 
1 


0 


which ensure zero direct and shear stress resultants and zero resultant moment along 
€=constant and =constant. Furthermore, the functions are orthogonal in the 
fundamental intervals (0 to 1) and therefore they can be used for expressing the Aity 
stress function F and also V*T as generalised Fourier expansions with the proviso 
that V*7, and its normal slopes, vanish at the edges. This restriction, however, 
can be disregarded for any practical calculations, as already discussed in Section |. 
The orthogonality relations for the characteristic functions “,(£) can be 


expressed as 
1 


=O when @ 
=1 when p=q 
where the symbol 4,,, denotes the Kronecker delta. 


It will now be assumed that the right hand side of equation (1) can be 
expressed as 


m=1 n=1 


Multiplying both sides of the equation (23) by “,, (€) VW, (y), integrating over the 
whole rectangle, and then using the orthogonality relations (21), 


bnn= ak | (é) (1)) dédy. (24) 


0 0 


Substituting equations (7) and (23) into equation (1), 


> { dun + = > Ann Pn’ =(), (25) 


m=l1n=1 b* a’ b* 
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Multiplying equation (25) by (£) and integrating, it follows that 


1 1 


0 0 


Equation (26) constitutes a system of linear simultaneous equations for the 
unknown coefficients a,,. For any practical problem it would be sufficient to take 
only the first few coefficients in the series (7) in order to obtain a reasonable 
accuracy. The resulting matrix formed by the linear equations (26) is particularly 
suitable for iteration, since the principal diagonal terms are larger than the 
remaining terms in any given row”. From equation (26) it is readily seen that the 
degree of convergence of the series in equation (7) depends approximately upon the 
factors (2,)~*. 


To facilitate setting up the simultaneous equations for the a,,, coefficients, the 
numerical values of the product integrals occurring in equations (26) have been 
calculated for values of m and p varying from | to 10 (see Table II). It can be 
shown that these integrals are given by the following expressions? : — 


1 


(B,* | P,P €=48,°B,," ¥mBm) [1 Pp (27) 


1 
0 


1 


When the a,,, coefficients have been determined, the thermal stresses are readily 
obtained from equations (8) to (10) with the aid of tables of the characteristic 
functions and their first and second derivatives. 


4, Expansion of V’7 in Terms of the ‘) Functions; ),,, Coefficients 
The b,,,, coefficients are evaluated from 


1 
ban =K | (2) V, (i) dédn . (24a) 
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The Laplacian of the temperature distribution V?7 can frequently be expressed as 
a product of two functions f and g, so that 


If f and g are given as a power series in ascending powers of £ and », namely 


f=Ac¢é ; g= 2 dak, ‘ (32) 


k=0 0 


With f and g expressed as Fourier series, 


cos cos krn 
sin sin 


it follows that 


\ sin | \ sin 


The integrals occurring in equations (33) and (35) have been evaluated for a range 
of values of A and m or n, and thus the b,,,, coefficients can be calculated with a 
minimum of effort. These integrals are given in Tables III, IV and V. 


5. Calculation of Thermal Stresses 


The practical application of the method developed in the present paper may be 
summarised as follows :— 


(i) Express V’T either as a power series or Fourier series. 


(ii) Calculate the b,,,, coefficients from equations (33) or (35), using Tables Ill 
or IV and V. 


(iii) Set up the simultaneous equations (26) for the unknown coefficients dns, 
using Table II. 


(iv) Solve the simultaneous equations by iteration (two or three iterations are 
usually sufficient for engineering purposes). 


(v) Calculate the thermal stresses from equations (8) to (10), using tabulated 
values for the <P functions and their derivatives. 
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TABLE VI 


THE @,,,, COEFFICIENTS OBTAINED BY SUCCESSIVE APPROXIMATIONS 


mn 


2nd iteration 


Coefficient Ist iteration 3rd iteration | 4th iteration 
10°a,,= 1°13630 1-15933 1-15997 1°15998 
,= 1:83064 1-88564 1:88683 1-88684 
354745 3°62817 3-62838 362838 
10°4,,= 2°36133 2:39913 2°39977 2-39977 
10°a,,= 8:50949 866991 8-67020 866921 
10°a,,= 9°44885 9°55214 9-55451 955456 
2°50808 2°48544 248429 248424 

= 3:25369 3:28309 3:28375 3:28376 


6. Numerical Example 


As a numerical example, consider a rectangular plate of aspect ratio a/b=2 
subjected to a parabolic temperature distribution in the y-direction, i.e. 


T (x, yy=4T. (y/b) 1 —(y/b)) 
where T, is the temperature at »=4. Assuming this temperature distribution, the 
practical steps outlined in Section 5 will be followed in detail. The Laplacian V*T 
is given simply by 
V°T (x, y)= —8T./b*. . ‘ 
Hence, from equations (32), 


f()=1 and g(y)= -87./b?=d, . (38) 


and the b,, coefficients are then calculated from equation (33), which in this case 
simplifies to 


0 


Evaluating the first eight coefficients 


b,,=0-690331 Kd,, by, =0-132328 Kd,. 

b,,=b,,=0-302242 by, =0-111356 Kd,, (40) 
b,,=0-192342 Kd,, = 0-084212 Kd,, | 
b,,=0-141051 Kd,, 
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It should be noticed that, due to symmetry of the temperature distribution, all 
coefficients for which m or n is an even number are equal to zero. Substituting 
a/b=2 in equation (26), the linear equations for determining the am» coefficients 
become 


1 1 


4 co 
0 0 


Hence, using Table II for the evaluation of the product integrals in the second 
term, equation (41) becomes, in matrix form, 


1215-05 —119-714 —93-6898 —119-714 94-6883 —75°1449 74-1042 —62-3316| a, | 
—119°714 4045-12 —299:567 94-6883 — 962-422 —282:766  236:944 —256:545]} 
—93°6898 —299:'567 153909 74:1042 236:944 — 467-871 —2568-91 —467:016] | a,, 
—119°714 94-6883 74-1042 30514-6 —962-422 59-4361 —753-203 49-3013] | a,, 
94-6883 —962:422 236944 —962:422 40844-6 223-655 —2408:32 202-915] | a,, 
—175:1449 —282:766 —467-871 59-4361 223°655 45777-4 370-064 —630-261| | a,, 
741042  236:944 — 2568-91 —753:203 —2408-32 370-064 664879 369388] | a,, 
—62:3316 —256°545 —467:016 49:3013 202:915 110401 a,, 

1-380662 

0:604484 

0:384686 
=| 0604484 

0:264656 

0-282102 

0:168424 

0:222712 


where, for the purpose of iteration, the equations are placed in such an order that 
the coefficients Gn, are decreasing. This has been achieved approximately by 
arranging these coefficients according to the ascending values of the principal 
diagonal terms. Since the principal diagonal terms are larger than the remaining 
terms, this matrix is particularly suitable for the application of iterative methods 
(see Ref. 7). The values of the @mn coefficients obtained in each successive 
iteration are given in Table VI, which shows that only two iterations would have 
been quite sufficient for engineering purposes. 


The thermal stresses can now be calculated from equations (8) to (10), using 
tabulated values for the characteristic functions and their derivatives. The pictorial 
representation of the distribution of the stresses o,, 7, and 7,, is shown in Figs. 2, 
3 and 4. A perusal of these distributions indicates clearly that the o, and ty 
stresses for plates of larger aspect ratios are confined only to the regions near 
the ends. 
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PicurE 2. Distribution of direct stresses +, due to temperature T=4 T(y/b) {1—(y/5)}; 
symmetrical about x=a/2 and y=b/2. 
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Ficure 3. Distribution of direct stresses ov, 
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due to temperature T=4 T(y/b) {1—(y/b)}; 
symmetrical about x=a/2 and y=5/2. 
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Figure 4. Distribution of shear stresses =, due to temperature T=4 T,(y/b) {1—(v/b)}; 


anti-symmetrical about x=a/2 and y=b/2. 


REFERENCES 

1. HELDENFELS, R. R. and Roperts, W. M. Experimental and Theoretical Determination of 
Thermal Stresses in Flat Plates. N.A.C.A. T.N. 2769, August 1952. 

2. SINGER, J., ANLIKER, M. and LEDERMAN, S. Thermal Stresses and Thermal Buckling. 
Wright Air Development Center Technical Report 57-69, Polytechnic Institute of Brooklyn, 
April 1957. 

3. Horvay, G. The End Problem of Rectangular Strips. Journal of Applied Mechanics, 
Vol. 20, pp. 87-94, March 1953 and pp. 576-582, December 1953. 

4. MENDELSON, A. and HirscHperG, M. Analysis of Elastic Thermal Stress in Thin Plate 
with Spanwise and Chordwise Variation of Temperature and Thickness. N.A.C.A. T.N. 
3778, November 1956. 

5. 


PRZEMIENIECKI, J. S. Temperature Stresses in Aircraft Structures. Ph.D. Thesis, Chapter 8, 
University of London, 1958. 


Youn, D. and Fetcar, R. P. Tables of Characteristic Functions Representing Normal 
Modes of Vibration of a Beam. Publication No. 4913, University of Texas, 1949. 


TIMOSHENKO, S. Theory of Plates and Shells. First Edition, p. 227. McGraw-Hill, 
New York, 1940. 


TIMOSHENKO, S. Vibration Problems in Engineering. pp. 331-336. Constable, London, 
1937. 


78 The Aeronautical Quarterly 


Tay / 
S 

b/2 
y 
| 

1 
hi 
sl 
tt 
a 
g 
il 
0 
h 
t! 
F 


of 


Supersonic Flow Past Slender Pointed Wings 
with “Similar” Cross Sections at Zero Lift 


W. T. LORD 
(Royal Aircraft Establishment, Farnborough) 
and 
G. G. BREBNER 
(Royal Aircraft Establishment, Bedford) 


SUMMARY: Some recent theoretical work on slender pointed wings at zero lift is 
co-ordinated and extended. The wings considered may have any pointed plan 
form shape, provided that the trailing edge is straight and unswept. The root 
section profile and cross-section shapes are arbitrary, provided that, on any 
one wing, the latter are “descriptively similar” (diamond or parabolic biconvex 
for instance), though not necessarily geometrically similar. The chief aim of 
the work is to find wings with simple geometry, low wave drag and pressure 
distributions which are unlikely to be seriously affected by viscous effects. 
Wave drag and pressure distributions are calculated by slender-wing theory. 
General formulae, which are both simple and instructive, are given for the wave 
drag and the overall pressure distribution, with particular emphasis on the root 
pressure distribution. Results for a number of wings of special interest are 
presented and discussed. 


1. Introduction 


Recent investigations into the properties of slender pointed wings have included 
some new work on the pressure distribution and wave drag at zero lift. Attention 
has been concentrated on isolated wings with sharp edges and simple cross-section 
shapes, unlike those which occur on wing-body combinations. The first paper in 
this new study, by Weber’, gives extensive information on the pressure distributions 
and local and total wave drags on wings of delta plan form of which all the cross 
sections are diamond-shaped. A method for designing thickness distributions to 
give low wave drag has been outlined by Lord and Green), but again the 
illustrations are confined to delta wings with diamond cross sections. The only 
other plan form so far considered is the so-called parabolic gothic plan form, which 
has a parabolic distribution of local span and a straight unswept trailing edge. 
Weber“ has discussed the properties at zero lift of such a wing with diamond cross 
sections. In the present paper this work is generalised to cover slender pointed 
wings with sharp edges, of arbitrary plan form and cross-section shape, provided 
that the leading edge is swept behind the Mach lines, the trailing edge is straight 
and unswept, and all the cross sections are “descriptively” similar, diamond or 
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parabolic biconvex, for instance. The primary purpose of this paper is to help the 
systematic development of wings with simple geometry, low inviscid wave drag, 
and theoretical pressure distributions which are not likely to be sensitive to the 
effects of viscosity, so that the inviscid wave drag may be substantially achieved, 
The design philosophy behind the investigation of slender wings is not discussed here 
(it has been treated in detail recently by Maskell and Weber’), nor is the problem 
of choosing a suitable basis for comparison between slender wings of different plan 
forms and thickness distributions. 


To study systematically the effect of wing geometry on wave drag and pressure 
distribution, a simple and reasonably accurate calculation method is required, so 
that a number of cases can be covered without undue labour. For this purpose the 
slender-wing theory has proved very valuable. It is a form of slender-body theory, 
with boundary conditions applied and velocities calculated in the chordal plane of 
the wing and not on the wing surface. This procedure is clearly questionable near 
the apex of the kind of wing considered here, but the region of greatest interest is 
near the trailing edge, where the thickness is small compared with the span and 
the theory is more likely to be valid. The pressure coefficient is taken to be 
proportional to the streamwise perturbation velocity, which is not a very close 
approximation at some parts of the wing. A guide to the effect of these approxi- 
mations at supersonic speeds may be obtained by considering, for example, slender 
bodies of geometrically similar elliptic cross sections®’. The results of the present 
paper for isolated wings show how the slender-wing theory suggests important 
qualitative trends and gives quantitative closed formulae which may serve as the 
basis for experiments. For the present, the effects of bodies and nacelles on the wing 
characteristics may be assessed from existing papers‘ * 


In Section 2, after transforming the co-ordinates to make the root chord, semi- 
span and volume equal to unity (the basis for comparison used throughout), the 
formulae for the pressure distribution and wave drag of wings with “‘similar”’ cross 
sections are derived. To appreciate more easily the effects of cross-sectional area 
distribution, Section 3 describes the geometry of wings with sextic polynomial area 
distributions, However, the achievement of low wave drag in viscous flow will 
probably depend on the pressure distribution, and in particular on the root pressure 
distribution, which may be satisfactorily adjusted by specifying root profiles rather 
than area distributions. Thus, in Section 4, a method of designing thickness 
distributions is suggested which involves examining only wave drags and root 
pressure distributions. The examples in Section 5 illustrate this method for delta 


and gothic wings. 


NOTATION 


G,,@2,@, coefficients in general sextic area distribution 
b,c, maximum span and root chord of wing, respectively 
f(x) non-dimensional root profile, defined such that fmax=f (1) = 1 


g(x) non-dimensional plan form, defined such that gmax=1 
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1 


h(y) non-dimensional cross section, defined such that | h(n) dn=1 


Amax maximum value of h(n), =h (0) 
j symbol for integer defining special family of root profiles 
m, position of root maximum thickness 
n position of maximum of cross-sectional area distribution 
q_ kinetic pressure 


1 
s(x) non-dimensional area distribution, defined such that s(x) dx=1 
0 


Smax maximum value of s(x), =s(n) 
t, root thickness/chord ratio of wing, =VZ, max/bc,” 
x non-dimensional chordwise co-ordinate, = X /c, 
y non-dimensional spanwise co-ordinate, = Y /4b 


z non-dimensional thickness-wise co-ordinate, =Z/{V /(2bc,)} 
Z, (x) non-dimensional root profile 
Zomax maximum value of Z, (x), (#1) =2, (x)/f (x) 
A, coefficient in general sextic area distribution 

B=(M?-1)'/? 

D_ wave drag 
F(x) function depending on s (x) occurring in ¢,’ 
G(x) function depending on g(x), =g” (x)/g (x) 
H(x) function depending on g(x), =2g’ (x)/g (x) 


1 


0 


I,,1,,1, constants depending on h(n). /,=T/; (0) 


1 


l 2 

J constant depending onh (7), = | (n,) h(y,) log 

1 2 
00 


M Mach number 
N=2n - (2 - 7n+7n’) 
S cross-sectional area 
T, coefficient of x” in polynomial form of s’” (x) 
U stream velocity magnitude 
volume 
X,Y  chordwise and spanwise co-ordinates, respectively 


Z_thickness-wise co-ordinate 
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5(y) function proportional to distribution of trailing edge angle 
¢ trailing edge angle 
7 normalised spanwise variable, = y/g (x) 
@ non-dimensional perturbation velocity potential, /(UV /c,*) 
x=4b/c, 
A wave drag coefficient, = Dc,*/(qV") 


perturbation velocity potential 


2. Wings with “Similar” Cross Sections 
2.1. GEOMETRY 


Consider rectangular co-ordinates X, Y and Z with the origin at the apex of 
a wing and with X as the chordwise co-ordinate and Y as the spanwise co-ordinate, 
Consider a slender wing which is symmetrical about the planes Y=0 and Z=0, 
and let the span, root chord, and volume be denoted by b, c, and V respectively. 
Let the area of a cross section of the wing be denoted by S and let < be the trailing 
edge angle of a chordwise section. Further, let the wing be in a supersonic stream 
in the direction of the X-axis, and let U, M and q be the speed, the Mach number 
and the kinetic pressure respectively in the undisturbed flow. Let “ be the 
perturbation potential of the flow past the slender wing, and let D be the wave drag. 


For convenience the co-ordinates are transformed to make the semi-span, root 
chord and volume equal to unity. The transformation is effected by writing 


Then, the plan form shape of the starboard half of the non-dimensional wing may 
be denoted by 
y=g (4), 1, 


where g(x) is such that g’ (x) is continuous for 0< x < 1, g (0)=0, and g (1)=guu 
=1. The thickness distribution of the upper starboard quarter of the non 
dimensional wing may be denoted by 


z=z(x,y), OSy<g(x), 


where z (x, y) is such that z(x, g(x))=0, z(1, y)=0. Then it follows that 


(x 0z 
|] zQ.ydy, s’Q)= dy, 
y constant 
° 0 


it was written while the author was at Imperial College, London. 


g(z) 


0 


0z 
(SZ); 
y constant 
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here s(0)=0, s’(0)=0, s(1)=0, and it is assumed that s’ (x) is continuous in 
0<xr<l. (The condition s(1)=0 is a consequence of the condition z(1, y)=0 
and is applicable here because only isolated wings with sharp trailing edges are 
being considered. For isolated wings with thick trailing edges and for wing-body 
combinations, the condition to be applied is s (1) 40, and the extension to this case 
is not difficult.) Also it follows that 


1 


If the non-dimensional root profile z=z (x, 0) is denoted by 
(x) 


and has the maximum value Z,max at x=m,, then the root thickness/chord ratio 
t, of the actual dimensional wing is 


VZeuse/be,*. 


By defining the variable » by 7= a 


and introducing a cross-section ‘shape function” h(n) which satisfies the conditions 
1 


h(l)=0, h(n) dy=1, then wings of “similar” cross sections may be defined in 


terms of the three functions / (n), g(x) and s (x) by 


s(x) 
2 () 


The first and second x-derivatives of the thickness distribution are 


02 s(x) s(x) (x) 
s” (x) 2s’ (x) (x) 


and, in particular, 5 (n)=s’ (1) h 
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If also h(n) is such that its maximum value on the wing occurs when »=0, 50 
that 4 (0)=/Amax, then the root profile z, (x) is given by 


Zo (2) = 


The root profile therefore depends on all three functions h(n), g(x) and s(x), and 
since s’(x) and g’(x) are continuous in 0<x< 1 it follows that z,’(x) is also 
continuous in 0<x<1. For comparison of different root profile shapes it is 
convenient to isolate the “‘shape function” f(x) associated with z, (x) by writing 


Zo max 


Then f(x) may be used as an independent function to replace s(x), in which case 


s(x)= 8 


and this is very helpful from the purely geometrical point of view and is considered 
later. However, the presentation of the aerodynamics is simplified if s(x) is 
retained as an independent function at this stage. 


2.2. OVERALL PRESSURE DISTRIBUTION* 


By slender-wing theory the non-dimensional potential ¢ on the plane z=0 
may be written as 


9 (x, ys 0)=9, (x, y) +9. (x), 
where 9, and ¢, are given by 


1 [a 1 
y)= yi | dy, + 58 (x) log x, 


6.0) = logaB- |, 


with x and B given by x=5/(2c,), B=(M? —1)'?. 


The pressure coefficient at a point on the plane z=0 is taken to be proportional 
to 9, (x, y, 0), where 


$2 (x,y, 0) (x, y) +420’ (x), 


*A more recent form for the overall pressure distribution is given in the Appendix. 
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l 0*z 


+ (x) logy, 


~— L[s” (x) log 3B - (0) log x - (x,) log (x - x,) dx, | 


The formula for (9,:)yconstant Can be written in terms of x, 7 instead of x, y, 
and, using the symmetry of z(x, y) about y=0, the result for a wing with similar 
cross sections is 


- 


where (x)= | h(n) log| dn, 


1 1 


0 


The formula for ¢,’ is not affected by the transformation from y to », of 
course, but an alteration may be made in it by writing 


z 


(x,) log (x - x,) dx, =(s” (x) - s” (0)) log x - F (x). 


0 


That this may be done can be shown easily in the case when s(x) is a polynomial, 
which is the case to be considered subsequently, and then, if 


N 
@= 


n=0 


Therefore 9,’ may be written 


(x)= [s” (x) log +F (x) |. 
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Thus it follows that if the functions G(x) and H (x), which depend only on 
g (x), are defined by 


then the overall pressure distribution may be obtained from the formula 


constant — { F (x) 
(x) +1, +9 (log 4xB + log? 


2.3. PRESSURE DISTRIBUTION 


Probably the most significant single pressure distribution is that along the root 
profile of the wing. This is obtained by putting »=0 in the formula just given. 


1 


Since | h(n) dn=1, it follows that J, (0)+J, (0)=-—1, and since A(1)=0 it 


follows that 2/, (0)+ 4/7, (0)+/,(0)=-—1. Therefore, the root pressure distribution 
is given, on writing J, (0)=J,, by | 


+ (log 3xB + 


2.4. WAVE DraG 
By slender-body theory the wave drag coefficient A is given by the expression 


| (x,) 8” (x,) log — | dx,dx, - (x,) log 


For wings of similar cross-section shapes the formula for A becomes 


A= { (x,) s” (x,) log | dx, dx, 2s’ (1) je (x,) log (, 


+97(1) log xB) } | 
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If h(n) is defined to be symmetrical about »=0 then 


1 
2 
J= (n,) A (m2) log dn, 


this quantity J defined specially for wings with “‘similar” cross sections is the same 
in this case as the quantity k defined by Lighthill“* for general slender wings with 
straight unswept trailing edges. 


2.5. DISCUSSION 


The overall and root pressure distributions and the wave drag depend essentially 
on the parameters x and B and the functions / (y), g(x) and s(x) and in this section 
some comments are made on the nature of the relations between them. 


The dependence on , and B is very clear and very simple. Only the product 
xB is involved, and only the function log ,B. Therefore, if the pressure distributions 
and wave drag for one wing at one Mach number, corresponding to a value (\B),, 
say, are known, then the results for a combination of wing and Mach number 
corresponding to some other value (\B),, say, are given, provided that A(y), g(x) 
and s(x) are unchanged, by 


($02) yconstant (Pax)y constant — (x) [log (xB), log (xB).], 
328” (2) {log (xB)» log 


Ay (1) [log (xB)> log (xB). 


It is noticeable that only in the overall pressure distribution is the effect of 
h(y) complicated. For fixed yB, g(x) and s(x), the effect of varying A(n) on the 
root pressures and the wave drag is very simple to calculate. Thus, if suffix a 
denotes the results for a particular cross-section shape, then the results for any 
other cross section, denoted by suffix b, are given by 


(Pv2)y=0 — (Par)y=o= (x) Loa) 


1 
Ay - (I) - J). 


For comparing different wing shapes the most important results are those for the 
wave drag and root pressure distribution, the latter being regarded as generally 
characteristic of the chord-wise pressure distributions elsewhere on the wing. These 
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results indicate whether the inviscid wave drag is low enough and whether it is 
likely to be realised in viscous flow. Since the effect of variations in A (1) on these 
results is easily determined and unlikely to be predominant, the choice of h (1) may 
be left until the later stages of the wing design, when overall pressure distributions, 
spanwise variation of local wave drag and lifting effects are being considered. The 
formulae for the wave drag and root pressure distribution are therefore examined to 
find out the effects of plan form, g(x), and area distribution, s(x), and whenever 
a precise wing has to be specified the cross sections are taken to be diamond-shaped. 
Results for two wings with parabolic biconvex cross sections are quoted in Section 5. 


For fixed xB, h(y) and s(x), the effects on the root pressure distribution and 
the wave drag of varying the plan form between cases denoted by suffixes a and }, 
say, are given by 


(90x) u=0 (Pax)y=0 = { [(Go (x) 4H,” (x)) (G, (x) 3H.’ (x))] 5 (x) 


+ [Hp (x) Ha 8’ +8” log 


A, A, 


While the expression for (¢1,)y<o — (@u»)y=) May not appear particularly simple, it is 
significant that it is independent of yB and h(y). The result for the wave drag is 
of special interest, since it predicts that the wave drag is completely independent 
of the plan form’shape. 


The effects of variations in area distribution s(x) on both the root pressure 
distribution and the wave drag are complicated, and no simple general rules can 
be given. Therefore, it is best to consider a definite family of area distributions. 
Polynomials are convenient for wings with sharp edges and, in the next section, 
the general family of sextic polynomials is considered, since this is expected to be 
adequate for present purposes. 


3. Wings with Sextic Area Distributions 
For isolated wings with sharp leading edges, for which s(x) satisfies 
s (0)=0, s’ (0)=0 and s(1)=0, it is possible to put 


S(x)=A,x? (L- x) +.a,x°), 


where, in order that { s(x) dx=1, it is necessary that 
0 


A,=12 (1+ a,+ 
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The derivatives of s(x) involved in the previous formulae are given by 
(x) =A, [x (2 3.x) (3 - 4x) + (4 - 5x) (5 - 6x)], 
s” (x)= A [2 (1 - 3x) + 6a,x (1 2x) + 4a,x? (3 - 5x) + 10a,x° (2 - 3x)], 
x” (x)=A,[ 6+ 6a, (1 4x) + 12a,x (2 - 5x) + 60a,x* (1 - 2x))]. 
The coefficients 7’, introduced in Section 2.2 are thus 


T,=6A,(a,—1), T,=24A,(a,-a,), T,=60A,(a,-a.), T;= - 120A,q. 
Hence F (x)= A,x { 6+ 6a, (1 -3x)+ (27 55x) + 2 a,x? (44 - 75x) 


The terms required for evaluating the wave drag of a wing with a sextic area 
distribution are 


(l)= - A, (1+a,+a,+4,), 


1 
Js (x,) log (, A, (5 + E Ut 


0 0 


4, Design of Thickness Distributions for Low Wave Drag 
4.1. GENERAL 


In Sections 2 and 3 the equations for the pressure distribution and wave drag 
of wings with “similar” cross-section shapes and sextic polynomial area distributions 
have been established, and in this section these are used to study the problem of 
achieving low wave drag. The first step in designing a wing to have low wave drag 
in viscous flow is to design one to have low wave drag in inviscid flow. No 
optimisation procedure exists, but slender-wing theory indicates that the inviscid 
wave drag of a slender wing depends predominantly on its area distribution, and 
it therefore seems reasonable to choose wings with area distributions which give 
low wave drag by this theory. However, the root pressure distribution of such a 
wing may indicate that adverse viscosity effects are likely to occur, and the procedure 
finally suggested involves specification of the root profile while aiming at a reasonable 
area distribution. The wave drag should vary only slowly with ,B if it is hoped 
lo achieve it in real flow. 
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4.2. WINGS WITH SPECIFIED AREA DISTRIBUTIONS 


The three most important quantities associated with an area distribution of 
given volume are the magnitude and position of the maximum area and the slope 
of the area distribution at the rear; in terms of s(x), the important quantities are 
Smax, M and s’(1). A low and slowly varying inviscid wave drag may be achieved 
by arranging that these parameters lie simultaneously in ranges roughly as follows: — 


1:65 <Smax<1°70, 0:50<n<0'55, 0<|s'(1)| <7. 


These values may be achieved by defining the coefficients A,, a,, az, a; in terms of 
Ny Smax and (1): 


A,N =420n' (1 - + —8n+7n?) Sax (10 - 38n+49n? - 21n’)s’ (1), 


a,A,N= 420n? (1 (2+n)- -2n-52n?+ 63n°) Smax — 
(1 (20 - 42n+35n')s (1), 


a,A,N =420n (1 (1 +2n) + (ul ~43n+7n? + 42n°) Smax + 


+n(1—n) (10—63n? + 70n’) s’ (1), 
a,A,N= - 420n(1-n)> - — 5n+5n?) Smax —7n (1 —n) (2—6n+5n?) s’ (1), 


N=2n(1-n)§ (2-7n+7n’). 


To assess whether viscosity is likely to have a marked adverse effect on a 
wing with low inviscid wave drag it is necessary at present to rely on a qualitative 
inspection of the inviscid pressure distribution. The viscous and compressibility 
phenomena most likely to increase the wave drag are boundary layer separation 
and shock waves forward of the trailing edge, and these are most probable when 
the inviscid pressure distribution contains large compressions. Attention will be 
confined to the root pressure distribution which is thought, in general, to be 
representative of other sections. The following example demonstrates an unfavour- 
able type of pressure distribution on a delta wing. 


It is well known that a parabolic body of revolution, pointed at both ends, 0 
that the condition s’ (1)=0 is satisfied, has a fairly low wave drag coefficient which 
is theoretically independent of Mach number, and so it is reasonable to investigate 
slender wings with the same area distribution as such a body, for they should all 
have the same wave drag irrespective of cross section and plan form shapes. The 
appropriate area distribution is 


s(x) =30x? (1 - x)’, 
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90 that A, =30, a, = - 1, a,=0, a,=0 and it follows that 


n=4, A=150/=. 


(The chief reasons for choosing this area distribution are that it is simple and that 
it corresponds to a very simple and much-investigated body of revolution, but it 
may be noticed that its value of smax is not within the range suggested for low wave 
drag. It is possible to design a pointed body with lower wave drag than this 
example, but for the purpose of illustrating wings with theoretically shockless 
compressions the present area distribution is satisfactory.) The cross sections are 
taken to be diamond-shaped so that h (y)=2 (1-7), A (0)=hmax=2, 1, = —3/2. The 
plan form is a delta, so that g(x)=x, G(x)=0, H(x)=2/x and then the formulae 
for the root pressure distribution and root profile are 


{ (4-73) +2 (1-6x+ 6x") log , 


(x) =60x (1 x). 


The root pressure distribution of this wing (designated wing (1)) is shown in 
Figs. 6(a)-6(c) for yB=0-1, 0-3 and 0-5. It is clear that for all values of xB 
there is a shockless compression over about the rear 0-4 of the root chord, the 
size of the compression decreasing as ,B increases. The flow around this wing has 
been discussed in detail by Weber”. 


The root profile of wing (1) is shown in Fig. 5 (p. 98). It has a cusped trailing 
edge, which is bad from the structural point of view. 


The results for wing (1) are typical of the behaviour to be expected for any 
slender wing with maximum span at a straight unswept trailing edge and with an 
area distribution with s’(1)=0 and n=4. The adverse pressure gradient at the 
toot and the inflection in the root profile both seem undesirable. Since they are 
interrelated, it may be that the elimination of the inflection in the root profile may 
lead to a root pressure distribution without a large adverse gradient. It is not easy 
to specify the area distribution so that the root profile does not contain an inflection. 
It is therefore suggested that the root profile should be the primary choice and that 
it should be adjusted, within the restrictions regarding inflections, to give a suitable 
area distribution. 


43. WINGS WITH PRESCRIBED ROOT PROFILES 


The prescribed root profiles should not contain inflections forward of the 
trailing edge. Also, the trailing edge angles should be as small as possible so that 
[x (1)]*, the coefficient of the log ,B term in the formula for the wave drag, is small 
enough to give plausible wave drag estimates over the range of Mach number for 
which slender-wing theory is valid. These requirements lead to the consideration 
of root profiles with the trailing edge as the point of inflection. A family of such 
profiles—the ‘‘j-family”—is now defined and investigated. 
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The function to be defined is the “shape function” f(x), introduced in 
Section 2.1. 


The “‘j-family”’ of root profiles is defined by 


(j+2)9t2G+n 


(1 x) (1 -(1 - 

where j is an integer, the maximum thickness occurring when x=m, defined 
by m,=1-—(j+2)-'/“+, and the trailing edge angle being governed by 


The cases j=0, 1, 2, 3, 00 are illustrated in Fig. 2 (p. 95). These profiles have 
sharp leading edges. The extreme case j=0, which is exceptional in the sense that the 
trailing edge is not a point of inflection, is the biconvex parabolic section and is 
therefore significant in its own right. The case j=3 is the highest value of j which 
can be taken if the combination of f(x) with a polynomial g(x) is to lead to an 
s(x) which contains no powers of x higher than the sixth, so that the sextic area 
distribution is adequate. However, the extreme case j=0O is of interest since it 
illustrates the trend of the profiles for higher values of j. 


The procedure to be followed in establishing a thickness distribution of low 
wave drag by the prescription of the root profile is to choose the profile with the 
lowest value of j for which the area distribution satisfies the given criteria for 
Ay Sax and s’ (1). 
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TABLE I (Continued) 


Note: The functions of / o(), 1,(), 1,(9) are as follows: 


Function Diamond Cross Section . Parabolic Cross Section — | 
1 | 
log (1—7?)'/2 + log "\+ log 

| 


—log (1—1?)!/2— tog (222). 


—72)1/2 1 
—log (1 —?)1/2 + 2 
+ — +n? 
1 1+ 
| 
| 0 
+ — +7°| 
| 


5. Examples 


To illustrate the effect of changes in wing geometry on the wave drag and root 
pressure distribution, some results are now given for nine different wings. The 
details of these wings are summarised in Table I and Figs. 1-4; the results are 
plotted in Figs. 5 and 6. Fig. 5 shows the variation of the wave drags with the 
slenderness parameter ,B and Fig. 6 the root pressure distribution for three values 
of xB. It must be emphasised that all the calculations are for wings with the same 
root chord, span and volume, and comparisons between the wings on some other 
basis might lead to different conclusions. 


Considering delta wings first, three of the examples have diamond cross sections 
and two have parabolic biconvex cross sections. The first three (wings (2), (3) and 
(4)) differ in root profile shape and, therefore, in cross-section area distribution. 
The root profiles are given by j=0, 1 and 3 respectively, and are successive stages 
in the attempt to reduce the wave drag while retaining a pressure distribution which 
might be substantially achieved in viscous flow. The simplest profile, j=0 (wing (2)), 
is the parabolic biconvex shape. This wing has been treated by Newby“®, and 
Weber’. Figure 6 shows that the pressure at the root becomes steadily more 
negative in going from leading edge to trailing edge. This is certainly likely to be 
realised in viscous flow, but the large suction on the rearward-facing surfaces leads 
to a high value of the wave drag (Fig. 5). By “‘flattening out” the root pressure 
distribution near the trailing edge it should be possible to reduce the wave drag 
without prejudice to the attainability of the theoretical distribution in viscous flow. 
A first step in this direction is to use the j=1 profile shape (wing (3)), and Fig. 5 
shows that this does reduce the wave drag. By going on to j=3 (wing (4)), the 
largest value consistent with a sextic area distribution on a delta wing, the wave 
drag is further reduced (Fig. 5), the corresponding pressure distribution being shown 
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TABLE I (Continued) 


Note: The functions of /,(n), are as follows : 


Function Diamond Cross Section Parabolic Cross Section _| 
1 
log (1—1?)'/2 +7 log (==) + log (3—n?) 
I 
(1 —n?)1/2 3 4 
2 


1 
|-log +9? log + — 

7 


0 


5. Examples 


To illustrate the effect of changes in wing geometry on the wave drag and root 
pressure distribution, some results are now given for nine different wings. The 
details of these wings are summarised in Table I and Figs. 1-4; the results are 
plotted in Figs. 5 and 6. Fig. 5 shows the variation of the wave drags with the 
slenderness parameter xB and Fig. 6 the root pressure distribution for three values 
of xB. It must be emphasised that all the calculations are for wings with the same 
root chord, span and volume, and comparisons between the wings on some other 
basis might lead to different conclusions. 


Considering delta wings first, three of the examples have diamond cross sections 
and two have parabolic biconvex cross sections. The first three (wings (2), (3) and 
(4)) differ in root profile shape and, therefore, in cross-section area distribution. 
The root profiles are given by j=0, 1 and 3 respectively, and are successive stages 
in the attempt to reduce the wave drag while retaining a pressure distribution which 
might be substantially achieved in viscous flow. The simplest profile, j=0 (wing (2)), 
is the parabolic biconvex shape. This wing has been treated by Newby“), and 
Weber’. Figure 6 shows that the pressure at the root becomes steadily more 
negative in going from leading edge to trailing edge. This is certainly likely to be 
realised in viscous flow, but the large suction on the rearward-facing surfaces leads 
to a high value of the wave drag (Fig. 5). By “flattening out” the root pressure 
distribution near the trailing edge it should be possible to reduce the wave drag 
without prejudice to the attainability of the theoretical distribution in viscous flow. 
A first step in this direction is to use the j=1 profile shape (wing (3)), and Fig. 5 
shows that this does reduce the wave drag. By going on to j=3 (wing (4)), the 
largest value consistent with a sextic area distribution on a delta wing, the wave 
drag is further reduced (Fig. 5), the corresponding pressure distribution being shown 
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20 


FIGURE 5. Wave drags A, wings (1)}—{9), 0°05 < xB <0°5. 


in Fig. 6. This has only a gentle compression near the trailing edge, which should 
be attainable in practice. 


The two delta wings with parabolic cross sections, wings (8) and (9), differ 
from wings (2) and (4) only in their cross-section shapes. From the root pressure 
distributions in Fig. 6 it is clear that, compared with wings (2) and (4), the decrease 
in pressure over the forward part of the section and the decrease in suction over 
the rear part should lead to smaller wave drags and this is borne out by the values 
of A for wings (8) and (9) in Fig. 5. The change of cross-section shape does not 
affect the calculated decrease of wave drag with increasing Mach number. This is 
very steep for wings (2) and (8) with j=0, and becomes progressively less steep as 
j is increased. 


The remaining three wings, (5), (6) and (7), have the parabolic gothic plan 
form with diamond cross sections, and differ only in their root profiles and 
consequently in their area distributions. The simplest root profile is the parabolic 
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biconvex shape (j=0) of wing (5) which has been studied by Weber“ and Peckham 
and Atkinson®"’, This wing has the same area distribution as wing (3), the delta 
wing with j=1 root profile. Consequently these two wings have the same wave 
drag, but their pressure distributions are different. The highest value of j consistent 
with a sextic area distribution and a parabolic gothic plan form is j=2 and this is 
chosen for wing (6). For small values of xB the pressure distribution has a greater 
compression near the trailing edge than any of the other wings considered in this 
section and a correspondingly lower wave drag. At higher values of xB this 
advantage is lost, since the pressure distribution does not change much with ,B, 
but the slow decrease of A with increasing xB is probably more realistic than the 
bigger decreases obtained with other wings. The remaining gothic plan form, 
wing (7), has the same area distribution as wing (4), the delta with j=3 root profile, 
and therefore the same wave drag. The root profile of wing (7) does not fit into 
the “7” category. The pressure distribution (Fig. 6) has no compression and it 
should therefore be possible to reduce the drag further if a small shockless 
compression can be allowed to occur at the rear of the wing. 


The relative magnitude of the wave drags in Fig. 5 should not be interpreted 
ds indicating the absolute merits of different wing shapes in achieving low wave 
drag. If a different set of geometric parameters were kept constant as a basis for 
comparison, the relative magnitudes of the wave drags might be different. However, 
the effects of varying any of the geometric parameters can easily be calculated from 
the results given in this paper. 
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Appendix 
It has been shown recently (August 1958) by J. H. Durran (Senior Mathematics 


Master, Sedbergh), while on a short visit to Farnborough, that the overall pressure 
distribution can be written in the following simpler and more enlightening form:— 


(92), constant (92), 29 + p (x, n) 


y (x)/8 (x) 
1 


{1, (n) - 1, (0)} = h(n,) log { | 7-4, |/,?} dn,. 


0 


Some details of the derivation of this formula and some new results deduced from it 
are given in R.A.E. Tech. Note Aero 2591. 
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